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EDITOR’S NOTE 


The papers presented herein have been derived primarily from transcripts taken at the Flight 
Mechanics/Estimation Theory Symposium held October 15 to 16, 1974, and April IS to 16, 
1975, at Goddard Space Flight Center. In order to achieve a uniform format, a considerable 
amount of editing was needed and in the interest of time, only a subset of all the papers 
presented were prepared for this publication. Abstracts of all the papers are included how- 
ever and attempts will be made to improve the handling of conference material so that in 
future proceedings all papers presented will be published. 
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N76-10168 


THE DETERMINATION OF ORBITS USING 
PICARD ITERATION 


Rajendra Prasad Mikkilineni and Terry Feagin 

University of Tennessee Space Institute 
Tullahoma, Tennessee 


The topic of this presentation is the determination of orbits using Picard iteration, v/hich is 
a direct extension of the classical method of Picard that has been used in finding approxi- 
mate solutions of nonlinear differential equations for a variety of problems. The application 
of the Picard method of successive approximations to the initial value and the two-point 
boundary value problems is given below. 

The initial value problem, 

z=F(z, t) 

z(t,) = c, n = l,2, ... 

is solved '^y means of the iterative scheme, 

z„(t,) = c. 

whereby the nth approximation to the state vector z is computed from the (n - 1) approxi- 
mation (beginning with some initial approximation, z^ft)). The computation involves per- 
forming a simple integration or quadrature at each stage of the iterative process. The con- 
stant of integration is chosen such that each iterate satisfies the initial condition, z (t, ) = c. 

Similarly, the two-point boundary value problem, 

z = F (z, t) 

g(z (t,),z(tj)) = 0, 

is solved by means of the iterative scheme, 

z„ = F(Vi-0 

g(z„(t,),z^(tj)) = 0, n = l,2,... 


/ 



where g denotes certain constraints or boundary conditions that must be satisfied by the 
solution at t, and tj . 

This paper presents an investigation of the suitability of this type of iterative scheme for tne 
problem of estimating orbits. In the estimation problem, we again have to solve the differ- 
ential equation for the state. However, we are now given a set of imprecise observations, 
y(tj ) j” , wfiich are, in general, nonlinear functions of the state. The problem is to find 
a z(t) which comes close (in some sense) to satisfying the observations. The estimation 
problem can thus be written as 

z = F (z, t) 

subject to the condition that z(t) minimizes some function Q, where 

Q = Q [z (t, ), z (tj) z (t„), y (t , ). }• (tj), . . . , y (t^)] . 

For instance, Q could be the sum of the squares of the residuals (the residual being the 
difference between the observed and ine computed values of the observations). Again, the 
iteration is set up, 

now choosing the constant of integration for z„(t) such that 

Q [z„ (t, ). (t j), (t„ ), y (t ,)...., y (t^ )1 

is minimized. As before, each stage of the iteration requires that a quadrature is performed. 
In order to perform the integration, we have chosen to approximate the right-hand side, 

F(Zj| , , t), by n series of Chebyshev polynomials, T^(t), and to integrate term-by-term. 

Thus, the approximation 

j=0 

is made where r = 1 - 2 - t, ) is normalized tme. The coefficients of the series are 

determined readily using the orthogonality properties from the relation 

k=0 

where tj = cos (jir/N). The integration can be performed by manipulating these coefficients, 
giving 

N , 
j=0 
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where 


aj=(bj., forj = l,2,...N 

and a,, = the arbitrary constant of integration that is determined from the requirement that 
Q be minimized. 

We will now consider the application of this technique to an orbit estimation prob:';;n., 
namely, the determination of the orbit of an earth satellite from observations o( the ange 
and range rate by several tracking stations. 

There are no external forces except the gravity field of the earth. The equations o' motion 
for the satellite are written as a set of six, first-order equations for the six state variables; 

X, y, z, u, V, and w. For instance, R,(t.) represents the range observed at t, from the fifth 
tracking station, and range rate is the time derivative of the range denoted by R,(t|): these 
are given. For the computer simulation presented here, we have taken the nominal values of 
the range and the range rate and added random noise. The standard deviations assumed are 
30 meters in range and about 50 cm per second in range rate. 

For this particular problem, the Q function should be considered. This is the sum of the 
weighted squares of the residuals; 

$ m 

-EE 

s=l 1=1 
*1 “l 

w., 

5=1 i=l 

where 

g,(tj) = |(x(t,)-x,(t.)l^ 

and 

x,(t,), yj(tj), Zj(t,) = the fifth station coordinatei a* t,, 

R,(tj) = the range from the particular tracking station, 

g, = the range computed using the piesent iterate, 

• 

Rj(tj) = observed range rate, and 

= weight for the particular observation. 
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We can see that is the function of the state at tj, and represents the coordinates of the 
station itself. The minimization of this function Q with respect to the integration constant 
is carried out using Newton's method. 

Table 1 provides the data for the three tracking stations we have assumed. The longitude 
and latitude, number of observations, and the total interval considered— about 1 500 sec- 
onds-are listed. The elliptic orbit considered has ar. ;centricity of 0.0557, with a semi- 
major axis of about 7178 km. This corresponds to a perigee at 400 km altitude, an apogee 
at 1200 km altitude, and an inclination of 20'* to the equatorial plane. 

Tab’s 1 

Station and Orbit Data 



Some of the results obtained are depicted in figure 1 , which shows the error as a function 
of time for different iterations. The initial guess is off by about 70 km from the true solu- 
tion, and with two iterations the error is brought down to something like 1 5 km. 

In order to see the convergence properly, the log of absolute error is plotted as a function 
of the iteration number in figure 2. The error is plotted for several points along the trajec- 
tory. The first point is T = 0, and the last point is T = 1 5 1 3. 1 . Jt is seen that the conver- 
gence is linear as was expected because of its relation to the classical Picard method. 

One of the disadvantages of this method is that the convergence is obtain;?d only for arcs 
of length less than ab it one-third of a revolution. Another problem is that of developing 
a procedure for con ng tv/o successive arcs. 


4 












TIME (>) 


Figure 1. Error as a function of time for different iterations 
(elliptic orbit). 


On the other hand, there are many advantages: For example, the method is quite simple and 
does not require the linear perturbation equations. The solution is in the form of a poly- 
nomial and is convenient to store. There is no interpolation required, should the solution 
be needed at some intermediate point. The error in representing this by a polynomial can, 
of course, be estimated by observing the last few terms of the polynomial. Also, the method 
is not sensitive to a poor initial guess. As we have seen here, the initial guess was off by 
about 70 km and even then the process will converge without difficulty. 

The question is, is this method reaUy more efficient than a more traditional approach to 
solving a problem? Actually, we have not compared it with the existing methods. However, 
we have some estimates on the time taken for this particular problem-it takes about 30 
seconds on the IBM 360/65. 


LOG OF ABSOLUTE ERROR 



Figure 2. Convergence of the method for an elliptic orbit. 
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A FINITE ELEMENT GRAVITY FIELD? 

J. L. Jiinkins 

University of Virginia 
Charlottesville, Virginia 


The paper details an approach for constructing a globally valid, piecewise continuous family 
of locall;' valid potential functions. The thesis put forth here is that the higher frequency 
terms of the geopotential can be more efficiently computed from such locally convergent 
approximations (typically three variable power series of order less than 5) than from any 
globally convergent gravity representation. This approach appears to be a step in the direc- 
tion of voiding the recent trend that the better model we have of the geopotential, the more 
expensive it is to integrate orbits with it! Numerical experiments conducted thus far confirm 
the validity of the approach and that acceleration errors of 0 (10''® km/s^ ) are achievable. 
The trade-ofis involved in selecting the finite element shape and size, and the order of the 
loca' approximations versus resulting accuracy, computational speed, and storage require- 
ment., are currently under study. 
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MEAN RATES OF THE ORBITAL ELEMENTS OF A SATELLITE 
PERTURBED BY A LENS SHAPED MASCON 


M. t. Ananda 

Jet Propulsion Laboratories 
Pasadena, California 


A set of mean orbital rates are computed for a satellite perturbed by a lens shaped mascon. 
A disturbing potential in terms of the orbital elements of the satellite and the mascon 
parameters is developed. The partial derivatives of the potential with respect to the orbital 
elements are formed. These partials are averaged over the period of the satellite orbit to 
eliminate the short periodic terms. The averaged partials are substituted into the variation 
of parameters equations to give the mean orbital rates. In the limiting case, when the radius 
of the lens shaped mascon reaches zero the mean orbital rate due to a point mass is ob- 
tained. The orbital rates developed by the method described here are compared against the 
rates obtained by numerical differencing. The method developed here is used to reduce the 
Apollo-15 and -16 subsatellite data for lunar farside gravity determination. 
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N76-10169 


ON THE FORMULATION OF THE GRAVITATI'^NAL 
POTENTIAL ns TERMS OF EQUINOCTIAL VARIABLES* 


Paul J. Cefola 

Computer Sciences Corporation 
Silver Spring, Maryland 


One year ago I presented a paper at this symposium that described some work we had oeen 
doing with the method of averaging. One of the aspects of that paper had to do with a par- 
ticular choice of nonsingular variables that were called equinoctial variables. In particular, 
we investigated analytical averaging techniques in the equinoctial coordinate frame. Some 
very good results were obtained, and we developed some expressions for fairly low-order 
perturbations. We developed explicit expressions for the disturbing potential for the zonal 
harmonics up to and the third-body harmonics up to (a/Rj , where a is the semimajor 
axis of the orbit and R, is the distance to the third body. 

After that we reconsidered the problem and decided that one of the shortcomings of that 
work was the lack of general expressions for the disturbing potential. This lack was a prob- 
lem for several reasons: First, serious problems were encountered in extending the results 
to higher degree terms. The need to extend the results occurs, for example, in mission anal- 
ysis, where the problem of including the effect of one h^er degree third-body term or in- 
cluding the effect of an additional zonal harmonic occurs frequently. Second, having explicit 
expressions for each term made software development complicated because, for each new 
.'•dditional software was required, and there was the possibility of errors at each sljp. 

These are some of the reasons why it was desirable to develop general results for the disturb- 
ing potential in this coordinate frame, and that is what 1 propose to discuss now. We are 
going to consider the gravitational potential and will use several special functions-the Legen- 
dre and associated Legendre polynomials and the , which are called the derived Legendre 
fun :tions. Recently there has been a fair amount of interest in these functions in the context 
of the gravitational potential. (For example. Pines, 1973.) The various disturbing potential 
expansions are listed below. 


* A significantly expanded version of this paper was presented at the AIAA Aerospace Sciences meeting in January 1975. 
Copies of the preprint are available from the author. 
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Disturbing Function Expansions 
Perturbation nth Term 

Third-body harmonics 


Zonal harmonics 


Arbitrary geopotential 
term (n, m) 


Here we are considering the third-body harmonics and the zonal harmonics, in particular. 

(It was originally planned to have some expressions for the tesseral case, but they are still 
in the process of being checked.) We have introduced the semimajor axis into these expres- 
sions, even though it is not required, because we are going to average these expressions and 
want to use existing results in two-body mechanics. For the third-body harmonics, the 4> 
angle will be the angle between the vector to the satellite and the vector to the third body. 

In the case of zonal harmonics, ^ will be the colatitude. 

In averaging these potentials, a two-step process is used. The purpose of step 1 is to obtain 
expressions for the Legendre functions in terms of equinoctial orbital elements. In step 2, 
the averaging is performed. 

One of the important questions in deriving potentials in terms of classical or equinoctial 
elements is how to get the fast variable motion into the potential. There are several different 
options, and they have a direct effect upon the amount of manipulation and on the com- 
plexity of the final results. Also related to these options is the level of accuracy that can be 
achieved. For third-body harmonics, there are really two options; We can use the mean 
anomaly for both the satellite and for the third body and end up with a very general potential 
This was done by Kaula in 1966. Or we can use direction cosines (relative to an orbital refer- 
ence frame) to get the location of the third body. If direction cosines are used, there is not 
quite the flexibility in modeling that is had with the Kaula approach, but the expressions are 
a lot more compact. In our work with the method of averaging, the use of the direction co- 
sines is quite appropriate. 

For zonal harmonics, we bring in the motion of the fast variable by using the orbital true 
longitude. For the arbitrary geopotential term, by which is meant tesserals, there are iwo 
options that are somewhat related to the two options that exist for third-body harmonics: 

We can use the mean variable for the satellite and Greenwich sidereal time, which is particu- 
larly appropriate for the study of resonance cases; or we can use the true longitude of the 



(C„mCOsmX + S^^sinmX) 
\ = a-d 
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satellite and Greenwich sidereal time, which appears to be most appropriate if we want to 
assume that the central body rotational angle is fixed during the averaging. Such a potential 
might be particularly appropriate for lunar satelUtes. 

Several expressions for cos ^ are listed below. These correspond to various attempts to 
develop the third*body potential. The first expression listed for cos ^ in terms of right 
ascension and declination of the satellite and the third body was used by Kaula. He then 
used the addition theorem and the definition of the inclination function to obtain expres- 
sions for the Legendre polynomials in terms of the classical elements. These results have 
great generality because we can impose all types of resonance constraints on them, yet they 
are also very complex. 

Expressions for Cosine ^ 


cos ill = a cos v + P sin V 

o = ^ P,P = ^ Q 


cos ^ = a cos u 0 sin u 

u = v + cj,o = ^ • N,P = ^ 'M 

cos ^ = a cos L + 0sinL 

L = v + t») + n,a = ^ = - g 

More recent work in third-body perturbations uses the concept of expressing third-body 
motion by direction cosines relative to an orbital coordinate frame. Kaufman in 1970 and 
Lorell and Liu in 1 97 1 , in particular, worked with the apsidal coordinate frame. Subsequently, 
Kozai worked with a coordinate frame where the orbital x-axis was along the nodal aossing 
and the third axis was in the orbit plane. Essentially, he used the argument of latitude as 
his fast variable. At about the same time, we used the true longitude as our fast variable and 
also used the direction cosines of the third body relative to the equinoctial coordinate frame 
(FGW), in which the x-axis points at the origin of the latitudes. 

One of the problems involved in working with Legendre polynomials and, in particular, with 
the direction cosines formulation, was that there did not appear to be an addition theorem 
for that case. In most of the work using that approach, the basic recuraon relationships 
have been used to generate the Legendre functions. There have not been general formulas 
for the required Legendre functions with the argument a cos L + j3 sin L. This is the same 
problem I reported on a year ago. Recently, we took another look at the addition theorem, 
looked at the expressions in terms of the direction cosines, and found that we could reformu- 
late the problem in such a way that the addition theorem still applied; the steps in that deriva- 
tion are shown in equations 1 and 2. 
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First, we introduce a phase angle, L', and deflne the trigonometric functions-cos L* and 
sin L'— to fit the form of the a«iUi>:ion theorem. Then a, 6 , otj, and 63 must be defined. 

a = L cos 5= 1 

03 = L' sin 6=0 ( 1 ) 

00563= \/a^ +0* sin 6 3 = -^ 1 - = 7 

Here those angles are treated as arbitrary quantities. By substituting equation I into the 
addition theorem; 

n 

PJcos^) = PJ0)PJy) + 2 ^ 

m= 1 


(n - m)! 


(n + m)! 


P..(0)Pnn.W“sn’(L-L'). 


( 2 ) 


In equation 2 , 7 is the direction cosine of the third body relative to the vector normal to the 
orbit plane. 

If we define some and polynomials in a and similar to the w • Pines did in 1973, 
the addition theorem can be simplified. The final result is 

n 

P„(cos = P (0) P„( 7 ) + 2 ^ P„J0) q,„(7) ic Jot./o cos m L + SJaJ) sin m Lj (3) 

m= 1 


where 


and 



(4) 

S.c.« = !.(«♦ jor 

(5) 


Because = 0 for n - m equal to an odd integer, gives even powers of the 7 function. 
It can be expressed quite well in terms of a polynomial in a and p. Again, it is finite, so we 
have a general form, a modified addition theorem, that is useful in this case. 

To give an idea of what these polynomials look like, the first one that appears in the third- 
body perturbation is evaluated as follows; 

3 

P-fo cos L + ^ sin L) = ( 37 ^- 1) +— ((a*- cos 2L + 2afi sin 2L] 

2 ' 4 4 

1 r 3 3 1 ( 6 ) 

= y y(o*+/ 3*)+— (o*-U*)cos2L+3o08in2L-ll. 

This checks against the results that were presented in our paper a year ago. It has been 
checked for several cases and gives exactly the results obtained previously through very 
long and arduous manipulations; in fact, some errors were found in those previous results 
through comparison with the general formula. 
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Anotlier major point of this paper is to emphasize the computation of the potential by 
using the recursion formula, and below are listed the recursion formulas (equations 7 to 10) 
tliat are appropriate to the computation of (o cos L + ^ sin L). 


(n+l)P„^,(7) 

= (2n + l)rP^(r)-nP„ ,(7) 

(7) 


(2n - 1) 7 Q„.,,„(7) - (m + n - 1) Q„.j _„(7) 

(8) 


= aCJa,fi).0SJaJ) 

(9) 


= P CJa.fi) 

(10) 


Equation 7 is the standard recursion for the Legendre polynomials, and equation 8 is the 
recursion for the derived Legendre functions (this is given in Ananda and Broucke, 1973, and 
Pines, 1*^73). Equations 9 and 10 are the recursions for the and polynomials, respec- 
tively. 


The final formula for the averaged potential also has the same and polynomials, with 
k and h as the argument, and the same recursions will still apply. We are interested in ob- 
taining averaged potentials, and we want to average with respect to a time-oriented variable, 
so we average with respect to the mean longitude. In equation 1 1 , we take the standard 
formula for the average of (r/a)" cos m v and recast it in terms of the true longitude: 


1 

lit 




e“"‘-dX 


H)”(^ 


(k + jh)"*F 


m -n - 1 


m - n 


2 


,m+ l,h^+k^ 


)• 


( 11 ) 


It should be noted that the formulas are a hypergeometric series and that the polynomial 
previously referred to as appearing in these integrals has the argument k + jh. The symbol 
j is the square root of -1 . 

We wanted to obtain a recursion formula for computing the right-hand sides of these integrals, 
and it happens that one has already been derived and is presented in Cook’s paper on the PROD 
program (1973). The formula for these integrals is as follows: 


i /'( t )'' (- t )' « 


+ ih)'"B|;;^j(h2+k2) (12) 


where 


Tift *»ft*l nic-1 »k (^+k)(8-k-l) + 1^2^ ftk + 1 

=l,Bg -BJ+ (h +k 

4ic (lc+ 1) 


13 



So along with (a cos L + 13 sin L), the hypogeometric series and the eccentricity poly- 
nomials can be reached recursively. 

The following is a result for the nth term in the third-body potential expressed in terms of 
equinoctial elements and the direction cosines: 



where 


n 

^ * k*) Qn. m (T') fC™ («.0) C™ (k, h) + («,/3) S„ (k, h)l 


m=l 



(n - m)! / n + m + 1 
(n + m)! \ m 


(13) 


All of the terms in equation 1 3 can be obtained recursively. In writing a program, we would 
probably initially compute the coefficients ^ and then store them as data at the beginning 
of the program. We do not obtain all of the terms indicated in the equation. We get approx- 
imately (n/2) terms in that summation because some of the associated Legendre functions 
with zero argument are zero, depending on whether n - m is an even or odd number. We have 
a finite, closed-form expression, and we verified several of the results presented a year ago 
with this formula. The next step is to show how the same analysis is done for zonal har- 
monics. The potential due to the nth zonal harmonic is 

II IRX l>f" 

(t) 

where rl/ is the colatitude. Cos ^ can be expressed in classical elements by 

cos - sin i sin (v a>) . (15) 

In equation 1 5, i is the inclination and v is the true anomaly. We can express cos ^ quite 
well in terms of the equinoctial elements as is shown in the next equation. 

-2p 2q 

cos\l»= cosL+ sin L (16) 

1 + p* + 1 + p* + q^ 

In equation 16, p and q are equinoctial elements. If we assume auxiliary variables defined 
by 
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-2p 


and 




(17) 


2q 

l+p"+q^ 


the argument of the Legendre function in equation 14 has the same form as it did in the 
third-body case, and we can use the addition theorem we derived previously for (o cos L + 
sin L). 


For averaging over the mean longitude, we have to somewhat modify our averaguig formula. 
We could use the previous averaging formula with negative values of n, but that would result 
in an infinite series in the eccentricity, which can be seen by examining the hypogeometric 
series involved. Essentially , we made a quadratic transformation in t he hypogeom etric series, 
changed the arguments a little, introduced a variable x defined as 1/y^l -h^ -k^ , and found 
that we could get a finite series in the eccentricity and x for the required integrals: 


lit 



x*" ‘ (k + ih)™ B™(h^ +k^). 


(18) 


We also use the same polynomials, the polynomials in this case, which can be used in 
computing the potential for both third-body and zonal effects. There will be another saving 
in computation if both effects are included. 


Next we give a general expression for the potential due to a zonal harmonic, that is, the 
averaged potential in terms of the same eccentricity functions, derived Legendre functions, 
and the same a-|3 polynomials; 


a® 


(0)P„(r)x"" ‘ (h^+k^) 


n n 1 ( 

£ *"• O'’ ♦ k’) Q„ „ C) (C„ («,«) C„ (k. h) t <«,» S„ <k. 

m=l 


where 


W. 


n, m 



(n - m)! / n - 1 
(n + m)! \ m 


Again, everything can be derived recursively; terms that are useful in the the third-body case 
are useful here, and terms that are useful for the lower order zonal harmonics can be used to 
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compute the higher order ones. Therefore, I suspect that a very efficient computation pro- 
gram can be built using this formula, with the added advantage of having only one formula 
rather than many to verify. 
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SOME RECENT INVESTIGATIONS IN NUMERICAL AVERAGING 

David Wexler 

Aerospace Corporation 
Los Angeles, Californki 


Our interest in averaging techniques steins from a requirement to conduct large-scale mission 
planning exercises on a system of computers that gets overloaded whenever si*ch a study 
arises. These studies have been conducted at Aerospace using a Cowell propagator simply 
because of the necessity for accuracy in a high drag environment. However, due to limited 
computer resources^ a project has been initiated to develop a better method, which will 
involve using an existing algorithm or developing a new one. 

This paper presents two topics, tne first of which concerns a feasibility study that was con- 
ducted using two programs, MAESTRO and PECOS, to, see if averaging could be accurate. 
The second topic concerns some continuing investigations and some new ideas, generally 
untested. It is assumed that the idea of averaging with respect to trajectory propagation is 
familiar to everyone so it wiL not be reviewed here. 

The Maestro program, developed by Chauncey Uphoff and Dave Lutzky of Analytical 
Mccimnics Associates, was obtained by Aerospace for a feasibility test. It was initially 
designed as an interplanetary mission analysis tool for the Radio Astronomy Explorer-2 
(RAE-2) lunar orbiter. Several choices of variables and propagating techniques are available 
in M\ESTRO. For instance, a Cowell integrator, an Encke integrator, and several variation- 
of-parametor techniques, deluding the numerical averaging of the Gaussian variation-of- 
panameter ( VOP) equations, are available. 

The forces in the program are given in the radial, circumferential, and polar directions. 
MAESTRO has a fairly simple exponential atmosphere and has a precision averaging startup, 
wiiich means that osculating initial conditions in several choices of reference frames can be 
input and the program will automatically provide the initial mean state. 

The Aerospace version of MAESTRO has been merged with a program called PECOS. 

PECOS, itself, is a test bed for orbit planning simulations. Tlic interplanetary calculations 
were stripped out of MAESTRO and the force evaluations from PECOS are now used by 
both MAESTRO and PECOS. 

The output can be plotted by an auxiliary pr a i available at Aerospace. In particular, 
overlays of PECOS p.^-ecision and MAESTRO averaging output can easily be accomplished for 
comparison purposes. Two versions of mean-to-osculatinir * ar'* u ed. 
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On output, the Kozai form is used, and within the averaging quar'rature, the Iszak form is 
used. There is no particular reason for this order except for the program’s historical develop- 
ment. Fortunately, Iszak is very good for use under the quadrature because it is nonsinguiar 
at low eccentricity. 

With present program limitations, only certain parameter sets and integration types from 
MAESTRO are available. For instance, the Cowell integration is not available. Only two out 
of the original eight techniques are being used. They include the integration of a parameter 
set which is nonsingular at low eccentricity. This set is similar to. but not exactly the same 
as, the equinoctial elements. The averaging technique uses the same olcaients. 

The initial conditions for the test case that were used for the feasibility study were chosen 
from ty^acal orbits of interest to Aerospace. They were initially compiled for Terry Harter 
to run at GSFC on a test using the Goddard Trajectory Determination System (GTDS). At 
Aerospace they were run by Stan Navickas in the Mission Analysis Department with the 
MAESTRO/PECOS software. 

Table 1 shows that the four cases here have perigees ranging from 67 to 92 n.m. and apogees 
ranging from 1 56 to 223 n.m. The inclinations are all between 90^ and 1 1 0.5°. The life- 
times range from about 4.S days to 2 months; case A-1 has a 2nnonth lifetime and case 
B-l has a 4.S-day lifetime, and those are the two ca.ses that will be discussed here. 

Figures 1 . 2, and 3 are time histories of elements for case A-l , which had approximately a 
2-month lifetime. The computer runs were limited to 7 days because it is so expensive to 
run the Cowell integration in PECOS. This is the orbit which does the least of the four 
test cases. 


Table 1 

Initial Conditions for Four Test Cases 


Case 

Perigee 

Hp(n.m.) 

Apogee 

inciinauon (deg) 

C„ A/W(fr /lb) 

A1 

92 

156 

96.43 

0.0007295 

(0.002099) 

\2 

85 

I6I 

94.55 

0.0007295 

(0.002099) 

B1 

61 

207 

96.57 

0.0062981 

(0.001286) 

B2 

69 

223 

110 5 

0.0062981 

(0.001286) 
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The semimajor axis is plotted in figure I . The oscillating line is the output from the Cowell 
technique in PECOS, and the line down the center is the mean element output from 
MAESTRO. There has been no correction to put the short-period variations back on. 

Figure 2 is the eccentricity from the same trajectory. Again the oscillating curve is the out- 
put from the Cowell integration translated into Kepler elements. Figure 3 is the ai^ument 
of perigee. Figures 4, 5, and 6 are for case B-1 , which dies ir* 4.5 days. The dark l*ne down 
the center is the mean output from MAESTRO. It can be seen thai the program tracks the 
elements all the way to the end of the lifetime. 

Figure 6 is the argument of perigee, and, as expected, the envelope of the oscillations grows 
rapidly as it approaches the end of the lifetime and the eccentricity drops. The orbit be- 
comes circular, so the argument of perigee begins to osdliate more and more wildly. But 
the mean argument of perigee is still being tracked and follows the center of the oscillation 
envelope. 

Figures 7, 8, and 9 give the comparison between the Cowell integration and the corrected 
average propagation at about 4 days into the case having a 4.5-day lifetime. The orbit is 
generally tracked. The phasing and the magnitude have tracked quite well at the end of 4 
days. I am very hesitant to interpret the small errors there, because we have a problem with 
the interpolator in PECOS/MAESTRO right now. The step sizes in the integration for the 
averaging process are 90 minutes, but the output is roughly every 2 minutes. 

if the output of mean elements were plotted on a finer scale, it would be seen tnat the result 
has a small wiggle. This wiggle is probably due to the interpolator. There is some suspicion 
that the errors in the osdliating elements are due to that wiggle, but that is still a qualitative 
judgment. The iinportant thing is that the orbit is being tracked. 

We int^erpreted these plots as proof that the averaging technique does work. At this point, 
there had been no attempt whatsoever to make this program efficient. MAESTRO was sim- 
ply merged into a larger system using some of the logic from that larger system, so timinp 
from these runs were totally meaningless. Using the results proving the feasibility of the 
technique, we decided to go ahead and attempt to modify MAESTRO to make it a stand- 
alone version that would be efficient, using essentially the same logic that it has always had, 
but using a more accurate force model than it had when it was developed for GSFC. 

The modifications and developments which will go into this further study will be discussed 
now. 

The purpose of the current investigation is to increase the accuracy a~-i, if possible, to de- 
crease the computation time; this will of course involve a compromise between the two. 
CuTently, we are developing MAESTRO as a stand-alone version with modification for 
accuracy in a high drag environment, to run on the CDC 7600. The next step will be to 
replace the independent variable with the true anomaly and to remove the fast variable. 
Another concurrent activity, being handled by Dave Lutzky of Vector Sciences, is treating 
the short period effects within the averaging quadrature and for output with a Fourier ex- 
pansion from which the coefficients of the harmonics are computed automatically and 
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numericaily. This is work which was developed under funds from GSFC. The results will be 
inserted into MAESTRO to see if they will hdp to make the program more efficient and 
accurate. 

Now I would like to discuss some of the work I am doing. It is all formal at this time, none 
of it has been implemented in the program yet. 

One of the disadvant^es of using time as the independent parameter is that the averaging 
period is always difficult to interpret. For instance, it is possible to use the mean motion, 
which is determined by the initial oscillating elements or by the initial mean elements, to 
determine the period. But, however it is done, it will probably be difficult to be consistent 
and to initially choose the right time. 

Another problem is that the equations for the averaged elements have a term in their exact 
form which depends upon the time rate of the period. This is generally referred to as the 
Leibnitz term, and it depends upon the time rate of the period, the values of the true 
elements at both ends of the averaging period, and on the mean element at the midpoint of 
the averaging interval. 

Still another problem arises if constant time steps are chosen; For instance, if the integration 
step is chosen as the initial period of the orbit, then, farther on down the line, the time step 
will no longer correspond to the period, and there will be drift with respect to perigee or any 
geometric event in the orbit. If the true anomaly is chosen as the independent parameter, 
none of the above disadvantages apply. Another effect is that we get automatic regulariza- 
tion of the orbit sampling. Points are automatically clustered near perigee where most of the 
action occurs for the high drag situation. 

It is also possible to make repeated use of terms in the mean-to-osculating transfcH'mations, 
particularly if step size is chosen properly within the averaging quadrature. Also, the mean- 
to-osculating formulas appear in closed form, at least to first order in . This is due simply 
to the fact that the disturbing potential can be expressed in a finite number of terms, de- 
pending upon the true anomaly. 

In changing the propagation algorithm, we wanted to choose a state vector which is nonsin- 
gular for low eccentricity and inclination. For convenience, we chose to use the same state 
that MAESTRO now uses, except for the fast variable. If the independent variable is now 
chosen as the true anomaly, then the angle-time relationship must be tracked through an 
equation other than the one we had before. If the new fast variable is itself chosen as the 
time, we then have a new differential equation; 

dt - 

— = h/r + perturbations, 
df 

Because this leaves us with large oscillations, it is a poor choice for averaging. Suppose 
that, instead, the parameter Q (Stem, 1960) is chosen; 


20 



Q = M + cj + S2-nt 
n = (-2c)^ III 


where e is the orbital energ>’. Then h determines the period of the true anomaly. The Q 
will always be a slow variable, and, in fact, for the Kepler problem, it is a constant. 

The equations for the new system are given by 


wherj L is any element 


IE dE dt 
df ~dt df 


(E,t,0 


A 

df 


E=F„(E.t,0 


The exact averaged equations are 


E = 


- r 

J 

f-W 


EfOdf* 


df 



dt 


Fg(E,r,t)df'~ 


To make averaging efficient, the force evaluation is accomplished by approximating the 
elements in the int^rand to the mean elements, plus a correction, as in the standard second- 
order averaging technique. Also, since the time no longer appears explicitly in the state, it 
must be the result of a similar correcting transformation. This means that the mean-to- 
osculating transformation is required at every step, as it is in all second-order averaging 
techniques. 

Formally, the mean-to-osculating transformations are simple to derive. The equations 
needed are indicated below in abbreviated form: 


E* = E + 6Ep 


If it is assumed that Af is 2rr/2K, then values of the tngonometric functions of the harmon- 
ics of the true anomaly can be tabulated and never recalculated, because the same values 
will be repeated in true anomaly from revolution to revolution. This should save a signifi- 
cant amount of calculation and make for an efficient mean-to-osculating transformation. 


In the equations, 


N 


dE 

— = CONST + 
df 



(A cos nf + sin nO 

' n n 
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and 


N 

6Ep =^^(-A^ sin nf + cos nO/n, 

I 

and are derived from the Lagrange planetary equations with the J ^ perturbation only. 
The ( J are slowly varying functions of the state vector and supposedly would not 
need to be recalculated more than once per revolution, and probably not that often. This 
depends upon the strength of the perturbations. 

The evaluation of dQ/df is fairly obvious, and it is used here as an example; The forcing 
function, , is (M + w + S2 - n) + (n - n) - nt. Time appears explicitly, which may create 
some difficulties, at least formally, although it is doubtful that numerically this will give 
any trouble. The factor dt/df, which appears in the integrand, is given by 

dt h 

— = — + perturbations, 
df 

However, since all the parameters in the state are slow (there is no fast variable), probably 
the two-body rate, or just h/r^ , is good enough for the function dt/df. That will probably 
require some experimentation, but, since there is no fast variable, there is no purpose in 
carrying higher order terms in this factor. The term h is related to the time derivative of 
the eneigy, and that relates (at least in its dominant term) to the time rate of the semimajm' 
axis due to drag for the cases where drag is the most impcvtant perturbation aside from . 
The quantity (n - n) is fH-oportional to the disturbing potential and (M cu -t- A - n) is 
given by the Gaussian VOP equations. The MAESTRO propagator will be used as the basic 
software tool in the implementation of these techniques. 

Still to be considered are techniques for calculating an accurate argument of perigee at low 
eccentricity. There are some problems because there is a mix of the true mean anomaly 
and the mean argument of perigee in one of the transformations; some experimentation is 
still required. 

The exponential instability in Q, which is due to the presence of ht, must also be investigated. 
This is a problem which will most likely be more apparent on computers with short word 
lengths. With the CDC computers and their essentially inflnite word lengths, we have tried 
experimenting to see when we can make the equations blow up in a simple form, and we 
have been unable to And the instability numerically. That may give us trouble, though, 
and it will require more experimentation to see if there really is a problem. 
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Figure 5. Comparison of PECOS and (MAESTRO eccentricity plots for case B-1 
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Figure 6. Comparison of PECOS and MAESTRO argument of perigee plots for case B*l. 
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Figure 8. Comparison of Cowell integration and short period recovery for eccentricity. 
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AUTONOMOUS NAVIGATION FOR ARTIFICIAL SATELLITES 


Ptanav S. Desai 

Computer Sciences Corporation 
Silver Spring. Maryland 


The term “autonomous navigation" refers to the possibility of providing a satelUte with the 
sufficient number and type of sensors, as well as computational hardware and software, to 
enable it to track itself. In other words, there is no ground tracking involved. 

There are two classes of such autonomous navigation: passive and active. Passive means 
that the satellite does not get cooperation from the ground or from other satellites; active 
means that the satellite does get active cooperation from either the ground or from another 
satelUte Uke the Trackir »ata Relay Satellite Systen RSS). This active coopera- 
tion could come in tl.e .adio »gnals, laser beams, o. jcher kinds of signals or beacons. 

The basic reason for usi- ^ .tonomous navigation is to reduce the necessity lor gr- und track- 
ing, thereby reducing the over, '-’ding of ground tracking faciUties and also reduemg the cost. 
There are also other technical reasons. For example, with a fast satelUte, if there is a gap 
in the ground tracking data set, especially if there are drag and other prominent effects 
present, an autonomous navigation system could increase the accuracy of prediction between 
the two data sets. Another reason could be that the reaction time for noticing changes in 
the satelUte orbit would be reduced by autonomous navigation, even if it is used as a 
backup to ground tracking. 

This work is a conceptuaUzation effort made by Don Novak, Paul Beaudet, and the author 
at Computer Sciences Corporation. The Uterature is not exhaustive, but it should be noted 
that Howard Garcia’s paper has a summary of sensors as well as some discussion of the new 
sensor interferometer landmark tracker. 

The following considerations are important in such a feasibility study: First of all, it is 
necessary to be aware of what types of sensors are available (or could be made available) on 
a satellite to help in autonomous navigation. Then the observability arising from combina- 
tions or configurations of these sensors should* be checked. In other words, it should be 
determined whether a given set of sensors is suffleient under various conditions for 
determination of attitude and orbit of the satellite. The accuracy of the selected system 
and its reliability should then be studied to determine that, should one of the components 
fail, the other components would be enough to back it up. The choice of sensors basically 
depends on the estimation algorithm, in that we might choose either a coupled attitude 
and orbit determination scheme or a decoupled scheme. The computational hardware is 
still another factor. 
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Some potential sensors for use in autonomous navigation are listed below, but this is 
by no means an exhaustive list. Some satellites that have used, or are presently using, these 
sensors are listed: 

• Inertial measurement unit— ATS-F, OAO-2, OAOC 

• Star mapper-ATS-F, CTS, OAO-2, OAOC, OSO-I, OSO-7, SAS-B, SASC, SSS-A 

• Magnetometer-AE, AEROS, GEOSC, OSO-I, OSO-7, SAS-A, SAS-B, SASC 

• Solar sensor-SE, AEROS, ATS-F, CTS, GEOSC, IMP-H, I, J, RAE-2, SAS-B, 

C, SSS-A, Nimbus 

• Horizon sensors-optical -IMP-H, I, J, RAE-2, ^S-A 

• Horizon sensors-infiared— AE, AEROS, ATS-F, CTS, SASC, TIROS, Nimbus 

• Interferometer landmark tracker— ATS-F 

• Scanner/camera-SMS, Nimbus, Landsat 

The first is an inertial measurement unit, which is a system of gyros and accelerometers 
for determining inertial attitude and inertial acceleration of the spacecraft; it has been 
used on the Orbiting Astronomical Observatory (OAO). The star mapper is probably the 
most accurate of the attitude sensors for determining inertial attitude and has been used 
on many satellites. The magnetometer determines attitude with respect to the magnetic 
field of the earth, or the central body, and if there is a good model available for the 
magnetic field of the central body, it indirectly determines the satellite attitude. The 
solar sensor comes in many varieties, but basically it provides angles to the sun from the 
spacecraft frame. Horizon sensors can be either optical or infrared. For example, an 
optical horizon sensor was used on the Radio Astronomy Explorer-2 (RAE-2) mission; but 
infrared is more common and is used in a large number of missions. The interferometer 
landmark tracker (ILT), a new type of sensor, will be discussed later. 

So far the discussion has been limited to attitude-type sensors; however, they are by no 
means the only sensors that could be used for autonomous navigation. We could consider 
using non-attitude- or non-nav^tional-type sensors, including meteorological cameras 
and scanners that could be used in a landmark determination scheme, in which there is 
currently an interest. 

There are two other sensors to be considered: the one-way Doppler and the image correlator 
(IC); however, the image correlator has not yet been put on board. The one-way Doppler 
would determine range rate to known radio stations on the ground or to a tracking satellite. 
The image correlator is an advanced version of the landmark determination-type scheme 
where a computer would determine, through pattern recognition, the direction cosines in 
the spacecraft axes to a known landmark. 
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It is also necessary to decide how to combine the sensors for autonomous navigation, and 
a decision has to be made as to whether attitude and orbit should be determined in a 
coupled Qr decoupled sch<‘me. There is an argument in favor of a decoupled scheme be- 
cause orbital parameters, on the whole, vary less rapidly than attitude parameters. There- 
fore, longer spans of data can be used for orbit determination than for attitude. A 
coupled scheme would not take into account this difference in the memory span require- 
ment of the two. However, an initial determination could very well be coupled, so perhaps 
the best method would be to first determine a coarse orbit and attitude in a coupled 
determination, followed by fine attitude determination, followed by fine orbit determina- 
tion. So, for this discussion, it is assumed that we are going to determine attitude first, then 
orbit. 

The ILT (figure la) has a pair of antennas (at A and B), and its functicn is to determine the 
phase difference between the received agnals that arrive at the antennas. The ILT must 
be initialized as to what radio signal frequency it is looking for, which means that the 
approximate position in (X'bit must be known. If the satellite is geosynchronous, orbit 
determination does not have to be done continuously. 



ANTENNAS AT A. B (BASELINE AB ' dl 
(•) 



Ibl 


Figure 1. Interferometer landmark trackers phase shift geom- 
etry and circle of ambiguity, (a) ILT phase shift geometry; 
(b| ILT circle of ambiguity. 
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We determine the path difTerence, d cos 0, by determining the phase difference, because 
the two are proportional, and this gives us the angle d between the direction line to the 
emitter and the baseline of the ILT. However, as seen in figure lb, there is a conical 
ambiguity left in the direction to the emitter, because that angle is all that is known. 

To reduce ambiguity, we also have another pair of antennas (not shown in figure) where 
the baseline is perpendicular to the first baseline; we will have an intersection of two 
cones and, therefore, will reduce the ambiguity to just two lines. A little further analysis, 
perhaps over a period of time or using multiple landmark determination, might reduce the 
twofold ambiguity as well. 

In figure 2, we begin examining some sensor configurations from the point of view of observ- 
ability. Suppose we were trying to determine the satellite inertial attitude using a star 
sensor that determines the direction to one star only, and then supplement that with a 
landmark determination scheme. In other words, the direction cosines in the spacecraft axis 
to one landmark on the central body are determined using either the ILT or some other 
means such as a scanner or a camera. 




Figure 2. Star sensor configuration, (a) One-star, single landmark 
symmetry; (b) poor resolution geometry. 
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In a scheme like this, the star determines the satellite attitude ambiguous to the ~oll 
around the axis from the satellite to the star. Normally, the roll ambiguity is further 
reduced by having a second star, but suppose that second star is replaced by a landmark. 

In that case, all that is known is the angle between the star direction and the landmark 
direction, and that does not give the absolute attitude of the spacecraft unless its position 
is known. Spacecraft position can be ambiguous to within that circle (figure 2a), in fact, 
to within the whole cone, which indirectly results in ambiguity in the attitude. That would 
be considered an untenable or unobservable condition. 

Suppose we tried to reduce the ambiguity by looking at the central body horizon and 
obtaining indirectly the direction to the center of the central body. This could remove the 
ambiguity, but occasionally there is a situation where the resolution is poor, because we 
are now on two circles that graze, and as can be seen in figure 2b, the graze is quite large. 

If the ILT station was not placed in the plane formed by the line to the center of the 
central body and the line to the star, then the circles would intersect at two points instead 
of grazing, and that would be a better resolution geometry. 

Based on this discussion, we can eliminate as untenable those configurations for attitude 
determination which use only one star and one landmark and such other combinations 
which are conceptually equivalent, for example: 

• One star— ILT/IC (single station) 

• Sun-lLT/IC (single station) 

• One star -central body horizon 

• Sun— central body horizon 

• Moon horizon-central body horizon 

• Moon ILT-centrai body ILT (single station) 

The following are weakly-orbit-coupled configurations for attitude determination: 

• One star— central body horizon -ILT/IC (single station) 

• Sun-central body horizon -ILT/IC (single station) 

• Moon horizon -central body horizon-ILT/IC (single station) 

• One star— ILT/IC (multiple stations) 

• Sun-ILT/IC (multiple stations) 

• Moon horizon-ILT/IC (multiple stations) 

• Moon ILT-central body ILT (multiple stations) 



We are still speaking of attitude determination, and we would like it to be orbit decoupled, 
but these configurations are weakly orbit coupled in the sense that the ILT. for example, 
would need to be initialized with an approximate knowledge of orbit so that we would 
know which station to tune to. This list includes a one star/central body horizon /one land- 
mark conflguration as well as a one-star /multiple station configuration for the same reason. 
In summary, to have orbit decoupled attitude determination, perhaps two star sensors 
would be needed, and they could have an associated inertial measurement unit to back 
them up. 

We will now examine the information that can be derived by determining the direction 
to the landmark in the spacecraft frame. Figure 3 shows the spacecraft position, the 
line from the spacecraft to a known landmark, and d and 7 , which are the latitude and 
longitude of the landmark. The unknown subpoint of the satellite is , 7 ^ ), the distance 
of the satellite is r, and the satellite-to-landmark range is p. The following equations are 
based on the fact that this line intersects this sphere, and they show that the direction 
cosines of that line ought to be 1 ,1 ,1 , in the earth frame of reference: 

X y z 


COS 0 cosy - pfi + r cos 6 ^ cos 7^ 

(1) 

cos d sin 7 = p£y + r cos sin y^ 

(2) 

sin d = pfij + r sin 

(3) 


p = -r C • r +-/(6 • r)^ -(1 - 1/r^) 


where 

8 • r = 8j^ cos 0 ^ cos 7^ + 6^ cos si” T'q sin 9 ^ . 

Equations 1 to 3 can be solved to get the range explicitly; therefore, for every landmark, 
the three equations reduce to two equations. If 6 and 7 , the position of the landmark 
on the earth, are known, then there are three unknowns: r, 9^ , and 7 ^ . If we are 
interested in orbit determination, and the attitude is already determined, then we have 
three unknowns and two equations per landmark, so that we need at least two landmarks 
for orbit or position determination. 

If we could determine the subpoint (6 ^ , 7 ^^) on a continuous basis, then, by studying a 
history of the subpoint, we could get a track of the satellite on the ground and its 
characteristics would suggest a period. By using Kepler’s Third Law, for example, the 
semimajor axis could be derived indirectly from the period. In that case, even one 
landmark is, in principle, sufficient for orbit determination. 
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Figure 3. Spaoecraft/landmark geometry. 


We have assumed that attitude was already known, but suppose we were trying to use 
this ILT for a coupled attitude and orbit determination, perhaps in a coarse way, for 
preliminary locking on. How many unknowns would then be had? 

In figure 3, we have these directional cosines to the landmark, which are supposed to be 
in the geographic frame. What we really would know from the satellite instruments would 
be the direction cosines in the satellite frame of reference. Therefore, indirectly there 
are three unknowns involved liere which represent the transformation from the satellite 
frame to the geographic frame, which essentially means the satellite attitude, that is, the 
three attitude angles of the satellite. We then have six unknowns— the three satellite 
attitude angles plus three positional unknowns— and two equations per landmark, so we 
still would need only three landmarks. If we had strategically located a sufficient number 
of landmarks, accounting for possible cloud cover, even then perhaps just a few would be 
enough for a coupled attitude and orbit determination. (Expressions for the sensitivity of 
this kind of determination Lave been developed in our report and are available for anyone 
interested in performing their error analysis.) 

In summary, it seems that by using a variety of sensors it is conceptually possible to have 
autonomous navigation. However, the details would have to be worked out for each type 
of orbit. For example, a geosynch onous orbit would require a different type of configura- 
tion, perhaps, than the 2-hour satellite, so we do not have a general conclusion or recom- 
mendation valid for all types of satellites. 

The second part of this paper presents a standardized autonomous navigation system 
(figure 4) which was designed by Computer Sciences Corporation personnel for possible 
future use in a standardized type of satellite. In the base are the computer, electronics, and 
some gyro hardware. On the two sides of that base are two gimbal units. The right4iand 
gimbal unit is a two-star sensor, located near the beam splitter, which would determine 
the satellite’s inertial attitude. The main navigational unit (lower left) consists of a gimbal 
unit with three arrays of infrared detectors arranged conically at the tip of the unit. One 
small section of that array is high resolution, and two outside sections are low resolution. 

This is supposed to determine three directions to three points on the horizon of the central body. 
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Figure 4. Possible configuration for the proposed standardized autorromous rravigation system. 

As seen in figure 5, three lines of sight to the central body horizon determine a uni ^ue cone, 
so that the satellite position is known. The following are the equations for Xj, y., < . , the 
central body coordinates in the spacecraft axes: 

h2,Xi + hjjy. + h„z, =(r*-R2)^ 


By using the information from the star sensors, the equations can be transformed into 
inertial coordinates. 

Figures 6 and 7 show the degree of accuracy we can get from this system. Depending on 
the angular accuracy of the IR detectors, different errors are obtained at different altitudes. 
At 1000 km altitude, and with a 0.1° angular precision, we have less than a lO-km error, 
in fact, close to a S-km error in altitude (figure 7). In figure 6, the horizontal component 
error, the curves are flatter, but again we have the same order of magnitude accuracy. 
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POSITION ERROR'^KILOMtTEHS 



Figure 5. Three lines of sight to the central body 
horizon determine a unique cone. 
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SPACECRAFT ALTITUDE^KILOMETERS 

Figure 7. Variations in altitude error 
as a function of altitude. 
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IDEAL ELEMENTS FOR PERTURBED KEPLERIAN MOTIONS 


A. Deprit 

University of Cincinnati 
Cincinnati, Ohio 


The motion is referred to Hansen’s ideal frame, its attitude being defined by its Eulerian 
parameters. The parameters selected to determine the motion in the orbital plane cause no 
singularities or indeterminacies for small eccentricities; they have been chosen with a view 
of making the right-hand members of the equations as simple to program as possible. 
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SOME RESULTS IN THE FUNDAMENTAL GEOMETRIC THEORY 
OF ONBOARD DIRECTIONAL SENSORS 


Bertrand T. Fang 

WolfPcscarcIi and Development Corporation 
Riverdale, Maryland 


Earlier in the symposium. Dr, Velez mentioned the need for standardizing the many differ- 
ent of attitude sensors. The primary attitude sensors now in use are directional sen- 
sors, which means that they sense or measure certain external reference directions relative 
to the spacecraft body. Interferometers, magnetometers, horizon scanners, and star track- 
ers are examples of such sensors, and this paper will take a very elementary look at some 
of their fundamental geometric properties, with the hope that, by studying the basic 
geometry, we can derive a set of equations which are applicable to different kinds of sen- 
sors. 

First we will look at the basic ingredients of directional information, the manner in which 
the observation equations govern these basic ingredients and convey attitude and orbit 
information, and also the manner in which errors enter into these equations. The second 
topic to be discussed is attitude observability; in other words, what combinations of these 
direction measurements will resolve attitude unambiguously. Lastly, we will consider some 
of the concepts developed to study horizon sensors. Horizon scanners are of interest be- 
cause they also present orbit information and are somewhat unique in that the measurement 
is not actually a particular vector, but rather a scanning vector, which is tangent to the 
spherical earth. 

In analyzing these different directional sensors very carefully, it is found that, although there 
are many different kinds of sensors, the directional information can be divided into two very 
simple types, the first of which is given by the scalar product of two vectors as shown in 
figure 1. 

For this measurement, we have the drection R, which is the radiation from the transmitting 
station, and also a spacebome antenna, which is a direction fixed onto the spacecraft. The 
measurement is the phase difference of the radiation arriving at the two ends of the antenna 
baseline and is given by the dot product of the reference direction R, which is a unit vector 
from the transmitting station to the spacecraft, and the unit vector K, which is the space- 
craft-fixed direction. This measurement may be represented as the following observation 
equation: 
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= cos Q ■ 


(i)(f)= 


R • K 


= (R'l |K‘| = IR‘| (A,^| [K»I 


The superscripts I and B refer to the inertial and spacecraft axes, respectively, and is 
the direction cosine matrix relating the spacecraft and inertial axes and containing all of the 
attitude information. On the left-hand side of the equation, y^ is the measured quantity: 
K® is a spacecraft-fixed unit vector and is therefore a known quantity. Therefore, the 
attitude matrix, , is the only unknown quantity in this equation. 




Figure 1. Attitude determination measurement as made by a short* 
baseline interferometer on the ATS*6 spacecraft. 


Figure 2 shows the second type of direction measurement, made by the digital solar aspect 
sensor used on the Atmospheric Explorer Spacecraft. The vector R represents the line- 
of-sight from the spacecraft to the sun; i, j, and k are orthogonal unit vectors along 
spacecraft-fixed directions; and solar aspect angle Eg is specified by the cross-product of 
R and K. This measurement may be represented by the following equations: 

Rxic R.J lR'riA,/Bl (J“I 

y^ = sin Eg = • i = = 

IRXKI n/i-(R*K)' n/i-UR'J^ (A,/b1 (K“j/ 

and 

(R'r(A,^i [i»i 

- cos Eg = 

>/i-((R'r(A,/Bi 
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Being scalar equations, this second class of measurement contains more information than 
the first class, although it is necessary that 

Yj +y* = I. 

It should be noted that and X3 are related in a more complicated nonlinear way to 
the attitude matrix, . 



Figure 2. Orbit determination measurement as made by the digital 
solar aspect sensor on the AE spacecraft. 


The various problems of interest may be classified as follows by referring to the meas- 
urement equations: In the usual attitude determination problem, only the attitude 
matrix, Aj^, is considered an unknown. In the orbit determination problem, only the 
reference vector, R, is assumed unknown. In coupled attitude/orbit determination, both 
of these quantities are unknown. It is known that there are three unknown parameters 
for the attitude and three unknown parameters for the orbital position; therefore, in 
theory, at least six equations are needed to determine orbit as well as attitude. One 
good point about this type of equation is that it is applicable in general to different 
kinds of sensors: whatever errors exist must enter into appropriate quantities in these 
equations. 

The first type of error to be considered is timing error: that is, instead of making the 
measurement at time t, the actual measurement is made at time (t + At). With the timing 
error considered, the fust measurement equation assumes the following form: 
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y, (t + AO = (R* (t ♦ At)| (t ♦ AOI |K®1 
s (R'(OI (A,^3(t + A01 1K“I 
- (R'(Ol (A, ^(01 ((K“l + At in“l X [K“|) 

It can be seen that, senerally, timing error does not have much effect on the reference 
direction sensed. The main effect is that, instead of measuring the component of R along 
the K direction, the measurement equation is measuring the component along a rotated 
K direction as a result of the spacecraft angular velocity, e^ . The following equation 
shows how errors in given quantities enter into the measurement equation; 

y, +Ay, = (lR‘| + lAR'l) (a‘^^1 (|Ki“ + lAKl“). 

On the left-hand side of the equation. Ayj could stand for instrument reading errors, 
biases, and so on. It is independent of particular instruments: that is, for different instru- 
ments, we may simply have a different bias, and so on. The second quantity. AR', 
represents the unco-tainty about the reference direction. There mew also be uncertainties 
because the instrument is not sensing the true direction; for instance, ionospho'ic refrac- 
tk>ns will result in terms like this. The last term. AK^ , represents the instrument align- 
ment errors. So it can be seen that equations of this sort will be applicable to all 
diffo'ent sensors, with any errors entering only into the d quantities. 

The next topic to be considered is attitude observability. The first basic principle we are 
concerned with here is that attitude is a relative notion, for although we generally refer 
to the attitude of spacecraft in inertial space, we could also refer to the attitude of the 
inertial space relative to the spacecraft. The information contained in one description is, 
of course, readily transformed to the other, but, conceptually, it often might be clearer 
to look at the problem one way rather than the other. 

Another basic fact is that an attitude has three degrees of freedom and requires at least 
three independent measurements for determination. In addition, if all three independent 
measurements are related to only a single reference direction, the attitude cannot be 
resolved without ambiguity. This should be obvious, since, for a rigid body, the position 
of any three noncolinear points (or equivalently, two noncolinear vectors) are re>^uired to 
determine its attitude. So if two reference directions in inertial space are known, the 
attitude of inertial space relative to the spacecraft, and, therefore, the attitude of the 
spacecraft itself, is known. For directional sensors, of course, attitude observability is 
equivalent to the observability of two reference directions. 

The following equation is an example of the analytical formulation of observability. The 
unit vectors e, , , and e^ represent three reference directions. 
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e, *e 3 = 7 

Assume that they are unit vectors and are known in inertial space. The components of 
these vectors along three spacecraft-fixed directions. ^ measured, that is, 

ej*r^ • Tj =n,C3 • f3 = p. 

Taken together, these represent nine nonlinear equations for the nine components of e^ . 

, and C 3 . The attitude observability is equivalent to the uniqueness of the real solutions 
of these equations, which are not very easy. Since we aie not concerned with error at 
this time, any measurement would correspond to an attitude. Thus there is no incon- 
sistency, and the only concern is whether there is no uniqueness of solution. 

In general, the analytical determination of observability is difficult. Figure 3a shows a 
graphical construction which can help a great deal in providing insight into these prob- 
lems. In the figure, vector e is a reference direction, and vector r is a spacecraft-fixed 
direction. A measurement of the component of e along the direction r will limit the 
reference direction to lie on a small circle, which is the intersection of a plane with the 
sphere. It would be very convenient to represent directions as well as attitudes on a unit 
sphere. 

Figure 3b shows the second type of measurement, mentioned before, which gives the 
angle 0. In terms of geometry on the unit sphere, it is a half great circle, which represents 
all possible directions with the same meridional angle 0. Since this is a half great circle, 
it does provide a little more information than the first type of measurement. 

Figure 4 shov/s the result of multiple measurements. It is assumed that two small circle 
measurements have been made, or that two components of a reference vector have been 
measured, which can be represented as the intersection of two small circles. In general, 
there is a twofold ambiguity, because after two components of a unit vector have been 
measured, the third component is determined up to an ambiguity in its sign. This is 
reflected by the two points of intersection as indicated in the figure. 

Figure S shows what happens with half great circle measurements. In general, the inter- 
section of two half great circles would completely define a direction. But these figures 
show that a combination of the second type of measurement, which is a half great circle, 
with the first type of measurement, which is a small circle, may still leave some ambiguity. 

When two reference directions are available, the results are as seen in figure 6. There are 
three small circle measurements which give two possible positions of the direction , as 
shown by the intersection of two circles. The reference direction must lie on a small 
circle centered at e^ as shown in the figure. In addition, A and B intersect, A and C 
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Figure 3. Graphical construction on a unit sphere, (a) Small circle representing all possible e which has 
component d along 7 ; (b> half great circle representing all possible directions with tne same meridional 
angle 0. 



Figure 4. Intersection of two small 
circles gives two possible directions for 
e; only components along linearly in- 
dependent directions are independent 
measurements. 

do not intersect, and Cj is known. Also, e^ has two possible positions, which means there 
is still some ambiguity in the attitude. 

The first conclusion arrived at through these arguments is that, generally, five independent 
small circle measurements are required to determine attitude. Three of the measurements 
define one vector completely, and the other measurement, together with the known angle 
between the two refeience directions, determines the other reference direction. Secondly, 
if we want to obtain attitudes from three measurements, then at least two of these meas- 
urements must be great circle measurements. A third conclusion is that, in general, three 
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Figure 5. Half great circle measurements, (a) The X^*r plane is perpendicular to the great circle plane, and 
e is determined uniquely; (b) two possible directions for e. 


three small circle 
measurements 


Figure 6. Three small circle measurements concerning two 
reference directions e^ and e^. 

half great circle measurements plus any additional measurement would result in attitude 
without ambiguity. 

It can also be concluded that attitude measurement based on three small circle measurements 
of three reference directions requires that the sensors not be on the same plane. This is easy 
to understand if we invert the role of the attitude of the spacecraft and the attitude of the 
inertial space; it can then be seen that this case is really the same as that covered by the first 
conclusion. 

Other conclusions are as follows; In general, the symmetrical placement of sensors in* 
creases the chance of attitude ambiguity. Although attitude observability is very difficult 
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to determine analytically, a mechanical model of a unit sphere can be constructed, which will 
readily resolve attitude observability without difficulty. Lastly, although sometimes more 
than three measurements are required to determine the three attitude parameters, four or 
five measurements available may also provide some redundancy for data smoothing. 

Figure 7 introduces the topic of horizon scanners. The reference direction for the horizon 
scanner is the local vertical, and when a horizon scanner measurement is made, it means 
that the scanning ray is tangent to the spherical earth. The information available is the angle 
between the scanning ray and the local vertical. The measurement equation is given by 


where 


cos 6 



= -R(r)-s(T) 


cos a sin 0 (r) - sin a cos ^ (r) cos j \Mr) + X (r) 


I 

')■ 


a = the half cone angle of the scanner, 

T = time of horizon crossing, 

^ = pitch angle of the scanning axle. 

= roll angle of the spacecraft, 

X = scanner roll angle, and 

h = spacecraft altitude. 

It is evident from the above equation and from what has been said previously that the 
horizon measurement does not convey any information concerning the yaw about local 
vertical. 


There i,. a meusar<.ment equation like the above for each horizon transit. Since the scanning 
spo:d is f:si. ii miy be assumed in a first approximation that neither the spacecraft altitude 
nor the • it Jce ras changed from horizon entry to horizon exit. It may then be deduced 
immedi itely fvc-n the two measurement equations at entry and exit that the spacecraft 
has a roll angle. 





and a pitch angle, 


0 = p +a 
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or 


0 = 7T - /i + a. 


where 


sin ^ = 



. . . . 

- sin‘ a sin‘ — 
a 


!4 


dA 

tau a = tau a cos 

a 


and AX is the earth width angle as seen by the scanner. 



S (r + At» 


S It), tcaiming ray 


The existence of two possible pitch angles is easy to understand if it is realized that the two 
horizon measurements are now nothing but two small circle measurements about the local 
vertical. The extension of this result to higher approximations with the consideration of 
small attitude changes is straightforward and. in that case, the horizon measurements will 
relate to the spacecraft attitude as well as attitude rate. 

Sometimes two scanners with different half cone angles a' and a" are mounted on the same 
axle. In that case, th'^ first approximation pitch angle becomes 
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0= tan 


-1 


. ... AX" 

sin a cos - sin a cos 

a a 

cos a* - cos a" 


The redundancy eliminates the pitch angle ambiguity and the need for altitude information 
and also provides smoothing for the roll information. 

Actual measurements are often transit times rather than scanning angles and earth width. 

The lime information is transformed to angular information using the scanning rate, and 
any bias in scanning rate will amplify with time. Therefore, without periodic reinitializations, 
the spacecraft roll, which is directly related to the scanning angle, cannot be determined 
accurately. On the other hand, the pitch angle is related to the earth width and is more 
susceptible to triggering errors. 
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DESCRIPTION OF THE ROTATIONAL MOTION OF A 
NONSYMMETRIC RIGID BODY IN TERMS OF 
EULER ANGLES, DIRECTION COSINES. 

AND EULER PARAMETERS 


H. S. Morton 

University of Virginia 
Charlottesville, Virginia 


The rotational motion of a nonsymmetric rigid body can be described by a set of three time- 
dependei.t Euler angles. In the torque-free case, the natural angles are the 1-2-1 or the 1-3-1 
Euler angles if the angular momentum axis coincides with the space 1-axis, the 2-3-2 or 
2-1-2 angles if it coincides with the space 2-axis, or the 3-1-3 or 3-2-3 angles if it coincides 
with the space 3-axis. In each case, the first angle satisfies a differential equation whose 
solution can be simply expressed in terms of theta functions, which can be readily computed 
with the aid of rapidly-converging series. One can then determine the nine direction cosines 
and/or four Euler parameters, whose transformation properties are more convenient than 
the Euler angles. The 2-3-2 or 2-1-2 angles offer certain advantages if the body 2-axis is the 
principal axis of intermediate inertia. The analytic torque-free solutions provide a good 
base for studying the motion of a non.svmmetric rigid body in the prese.nce of perturbing 
torques. 


54 



ATTITUDE CAPTURE PROCEDURES FOR GEOS-C 
G. M. Lemer 

Computer Sciences Corporation 
Silver Spring, Maryland 


The scientific objective of the GEOS-C mission is to perform experiments to apply geodetic 
satellite techniques to solid-earth physics and oceanography. A spaceborne radar altimeter 
will map the topography of the ocean surface with a relative accuracy of ±1 meter. 

To meet the objectives of the radar altimeter experiment, GEOS-C will be gravity-gradient 
stabilized in a nearly circular orbit at an altitude of 843 km. The orbital inclination will be 
1 15° to maximize experimental coverage in the North Atlantic Ocean. A 21.5-ft, extendable 
scissors-type boom with a 100-lb end mass will provide passive pitch and roll stabilization 
and a momentum which will provide yaw stabilization. An eddy-current damper and mag- 
netic coil are also provided. 

Because no active attitude control hard A^are is provided and damper time constants are 
large, a detailed capture strategy has been developed at Computer Sciences Corporation 
(CSC) in coordination with AVCO and the Johns Hopkins Applied Physics Laboratory. 

This strategy has evolved from many detailed simulations and requires real-time attitude 
determination support for the initiation of critical commands. 
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POTENTIAL GEOSYNCHRONOUS ORBIT/ ATTITUDE RESOLUTION 
USING LANDMARK DATA 

C. C. Goad* 

Wolf Research and Development Corporation 
Riverdale, Maryland 


The information content in data other than conventional radar tracking is gaining increased 
popularity. This paper presents a ariance analysis of reducing the orbit and attitude state 
from iandmarks (ground-truth) data extracted from images taken at geosynchronous altitude. 

It is shown that comparable accuracy can be achieved with landmark data alone when com- 
pared with current optimistic reductions using cov' cntional types of data. 


‘Currently at National Oceanographic and Atmospheric Administration, Rockville, Maryland. 
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DETERMINATION OF ORBrTAL POSITION FROM EARTH 
AND SUN SENSOR DATA ON SMS-A AND IMP-J 


H. L. Hooper and M. A. Shear 

Computer Sciences Corporation 
Silver Spring. Maryland 


Attitude data from earth and sun sensors on SMS-A and IMP-J were processed in an attempt 
to refine the orbital elements as determined from tracking data. The results were checked 
a^inst additional tracking data and the discrepancies were investigated. The results have 
implications for the accuracy obtainable with any orbit-dependent attitude sensor. 
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ON-LINE ORBIT DETERMINATION AND ESTIMATION 
FOR ATS^ FROM REDUNDANT ATTITUDE SENSORS 


T. S. Eni^r, Jr. 

Btjsmessand Technology Systems. Inc. 
Seabrook. Maryland 


ATS^ is equipped with an onboard, two^xis interferometer which can provide direction 
cosines to earth-based transmitters. In addition, the spacecraft carries an earth scanner which 
can be thought of as an additional two-axis interferometer with transmitter at the earth's 
center. From the six-direction cosines thus available, both position and attitude determina- 
tion can be performed. This paper describes a procedure proposed for use both in the 
SAPPSAC experiment on ATS-6 and in the ATS-6 on-line attitude determination program 
at GSFC v/hich decouples position determination from spacecraft attitude. The resulting 
position pseudo-measurement is used in a constant gain Kalman filter fur estimation of orbit 
state. Propagation of the state estimate is accomplished with circular orbit perturbation 
equations. 
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THE OPTICAL SUT SENSOR AS A STANDARD SENSOR FOR 
SPACECRAFT ATTITUDE DETERMINATION 


James Wertz 

Computer Sciences Corporation 
Silver Spring, Maryland 


The idea for using an optical slit sensor as a standard sensor for spacecraft attitude deter- 
mination arose during a brainstorming session of the First Goddard Standardization Meeting 
on June 6, 1975. This paper describes the basic concept of the slit sensor, indicates what 
information is available from a single sensor and from two sensors, describes one possible 
standard sensor package, and compares the standard sensor package with the attitude 
package flown on the first Synchrrmous Meteorological Satellite (SMS). 

The slit sensor concept is one which has exciting potential as a standard attitude sensor 
for an enormous variety of missions, specifically, any mission using a spinning spacecraft 
or where rotating sensors or mirrors could be used. At present we are still in the stage of 
analytic studies-no experiments or design studies have been made. However, past exper- 
ience suggests that such sensors are feasible and should be relatively easy and cheap to 
build. 

The basic idea of a slit sensor is simple, as shown in figure 1 . It consists of one, or perhaps 
two, narrow slits with a 180° field of view capable of triggering on both the earth and the 
sun and distinguishing between them. There is no angle measurement per se. ilrere is simply 
a voltage change or pulse as the sensor crosses the sun or enters or leaves the d: ° the 

earth. The slit scans the sky by being mounted parallel to the spin axis on a sp... ig space- 
craft, by rotating itself, or by looking into rotating mirrors on a nonspinning spacecraft. 

The major advantage of such a sensor is that it sees the entire sky. If the earth and the sun 
are visible, it will see them. Equally important, however, is that the slit sensor returns a 
wealth of attitude data that is both easy to interpret and particularly amenable to sophis- 
ticated analysis techniques. 

It is important to keep in mind that the idea presented here is a standard sensor concept 
rather than a single piece of hardware; that is, we are interested in standardization of style 
and in manner of use. This implies great cost reduction, with little or no loss, and, in fact, 
possibly some gain in the versatility for individual missions. 

The greatest economy probably comes from the fact that there is only one ground proces- 
sing package, thus reducing the development costs. On the other hand, the experience 
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base that is gained provides economy, reliability, and confidence. There will be some degree 
of package duplication, and this will certainly reduce hardware costs for some missions. At 
the same time, the mission planner is free to adjust a variety of physical parameters to meet 
his needs, or to take advantages of advances in hardware or software design. In addition, 
the mission attitude analyst is free to concentrate his efforts on improving quality and 
accuracy of results without constantly having to start over by developing new models or 
combining a variety of old models for each new mission. 
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Figure 1. Optical slit sensor. 

It has been indicated that a wealth of attitude information is available from the slit sensor 
or sensor package. Next I will briefly describe the type of information obtained, and then 
compare it with that available from sensors actually flown on the recently launched SMS-A. 

In examining the information available from a single sensor as shown in figure 2, the most 
obvious data is a dihedral angle from the sun to the earth, measured to the earth midscan 
(D). However, this could also be measured to earth-in or earth-out for purposes of redun- 
dancy. A second type of data is the nadir angle, or the angle from the spin axis of the space- 
craft to the center of the earth, which is available from earth-width measurements. Figure 3 
is a plot of nadir angle versus earth-width angle, with nadir angle along the vertical axis 
and earth width along the horizontal axis. As can be seen from the figure, as the earth 
moves toward the sensor poles, it will subtend a larger dihedral angle. Therefore, the eaith- 
width angle can be used as a measure of the nadir angle. 

Tliere are two regions where there is some difficulty with this measurement: In the vicinity 
of the spin plane, the earth-width measurement is relatively insensitive to the nadir angle, 
and, in the vicinity of the poles, the sensor never leaves the disk of the earth. Thus, there is 
not particul. y good attitude information in two regions-along the Equator an& in the 
region of the poles. However, there is a large region in between where high quality attitude 
information is available. 
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Figure 2. Irtformation available from a single sensor. 

Both of these measurements are available from just the presence of simple pulse. If the 
total illumination following on the slit sensor is monitored, then additional information is 
available which can provide the nadir angle in the vicinity of the poles. That is. specifically, 
if the total illumination output from the slit sensor is monitored, then this will be a sinus- 
oidal oscillation. The amplitude of the sine is a measure of the nadir angle, and the phase 
of the sine relative to the sun is a measure of the dihedral angle from the sun to the center 
of the earth. 

This information alone is sufficient for attitude determination; therefore, we could stop 
with a single sensor. However, the information would be of relatively poor quality at nadir 
angles near 90“, which are relatively common, and near the poles .he information depends 
on intensity measurements. Bias determination and in-fli^t calibration would both be 
difficult. 

Thus, it is worthwhile to consider the nature of information available from a second sensor. 
In particular, if the second sensor is mounted at an angle to the first, rather than parallel 
to the spin axis, then there will be a great deal of additional information rather than just 
redundancy. Slightly different information would be available, depending on whether the 
second sensor looks into the same hemisphere as the first or is canted somewhat so it can see 
slightly up or down. For the sake of concreteness, the results of a 360°sensor will be con- 
sidered, but this is by no means necessary to take advantage of the second sensor. 
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Figure 3. Nadir angle versus earth-width dihedral angle for 
vertical slit sensor and 70^ diameter earth. 

There are measurements available from a second sensor that are similar to those from the 
first sensor, but there are important differences as well. The dihedral angle from the sun to 
the earth is measured between two planes that do not contain the spin axis, and the nadir 
angle from earth widths (shown in figure 4) comes closer to the poles. 

Specifically, if p is the angular radius of the earth and Q is the angle at which the second 
sensor is tilted to the first, then the nadir angle measurement will come within of both 
poles. Thus, full sxy coverage in the nadir angle is available by tilting the second sensor 
at an angle equal to the angular radius of the earth at whatever distance the spacecraft is 
operating. The intensity measurements are the same at the pole as well, if we wish to use 
them; and the same problem exists around the Equator where the earth width is relatively 
insensitive to the nadir angle. 

However, in addition to these two measurements, there are two measurements that are 
available from two slit sensors which are not available from either sensor singly. In par- 
ticular, there is a sun angle from the sensor crossing times, and, as illustrated in figure 5, 
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Figure 4. Nadir angle versus earth-wirlth dihedral angle for 
360° field of view slit sensor 35° from vertical and 70* 
diameter earth. 

there is a third independent measure of the nadir angle from the midscan crossing times. 

Here again the nadir angle is the vertical coordinate and the midscan-to-midscan dihedral angle 
is the horizontal coordinate. It can be seen in the figure that there is now particularly good 
data in the vicinity of the Equator, with good distinction of the nadir angle. In addition, 
since the two hemispheres are different, we no longer have the problem of ambiquity about 
two possible solutions. Therefore, one needs no a priori information at all in order to use 
this procedure to find the attitude. 

The following is a summary of the measurements that are available from two nonparallel 
slit sensois: There is one sun angle measurement, two independent measurements of the 
sun-to-earth midscan dihedral angle, and two earth-in measurements and two earth-out 
angles available lor redundancy purposes. There are three independent mea. urements of 
the nadir angle. In addition, if we are willing to monitor the total illumination falling on 
the slit sensor, then there will be two additional independent measurements of the nadir 
angle available. So it is seen that two slit sensors provide essentially full sky coverage 
with abundant attitude data. 
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Figure 5. Nadir angle versus dihedral angle from midscan of 
vertical slit sensor to midscan of 35° off vertical slit sensor 
for 70° diameter earth. 

Given the general concept of attitude determination with a pair of optical slit sensors, it 
should be possible to incorporate these sensors into a standard attitude package. The key 
to a standard sensor package is versatility; that is, a variety of sensitivity levels, a variety of 
spectral regions, and a variety of angular orientations. For example, the most likely spectral 
region for normal operation would be the infrared, since this provides well-defined earth 
horizons and intensity levels such that a single sensor could trigger on both the sun and the 
earth. At the same time, it is desired to incorporate the possibility of other intensity levels 
and other spectral regions for triggering on different celestial objects. 

There are many possible configurations for a standard sensor package. One possible configu- 
ration, shown in figure 6, consists of three sensors mounted on a single plate: One sensor 
would be parallel to the spin axis of the spacecraft and two would be tilted at adjustable 
angles. In normal operation, the vertical sensor and one of the others would be used for 
attitude, and the third sensor used purely for purposes of redundancy. Any two sensors 
could provide high quality attitude information over all, or nearly all, of the celestial sphere, 
and any one could provide adequate information over most of the celestial sphere. 
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Figure 6. Slit sensor package. 

So far we have discussed the general characteristics of slit sensor behavior Now we will 
compare the analytic performance of a slit sensor package with the attitude package flown 
on the first SMS, which was launched on May 19, 1975, into a transfer orbit and was 
eventually placed in a circular synchronous orbit near the Equator. It should also be 
pointed out that, relative to Mr Goad’s presentation, we are interested here only in the 
attitude sensors themselves, and not attitude determination from the visible and infrared 
spin scan radiometer (VISSR). 

There were seven attitude sensors actually flown on the SMS: two sun sensors, each with a 
field of view of 120^and five earth sensors: two primary earth sensors at 4 above and below 
the spin plane that were used for attitude determination in mission orbit and three other sen- 
sors that were used primarily for attitude in the transfer orbit. We will compare this with a 
slit sensor package consisting of three attitude sensors: one sensor that is par.^Ilel to the spin 
axis and two which are tilted, one at 8.5° and the other at 30° in the opposite direction. 

Accurate spacecraft attitudes require that sensor biases be accurately known. Therefore, it 
is of interest to determine how many biases would need to be measured to use the sensor 
information with precision. In the package flown there are 20 total biases (two plane tilt, 
six azimuth, five triggering level, and seven elevation angle), which, as a purely practical 
matter, tells us that some selection will be necessary. There is essentially no chance of 
determining 20 biases from the data that are available from the spacecraft. The slit sen ior 
package has a total of eight biases (three plane tilt, two azimuth, and three triggering level) 
that would need to be determined for all of the sensors to be utilized. 


65 



In normal circumstances, however, all the sensors in either package will not be fully utilized, 
and all of the biases will not need to be determined. Specifically, in both cases, a minimum 
of five biases is required for accurate attitude determination. However, it should also be 
pointed out that these five minimum biases for the sensor packages flown provide accurate 
attitudes only near the mission orbit and, in particular, do not provide accurate attitudes 
during most of the transfer orbit or during any other maneuvers or mishaps that might occur. 
On the other hand, since one sensor is fully redundant in the slit sensor package, the five 
minimum biases there provide accurate attitudes over the entire celestial sphere, including 
both the transfer orbit and the mission orbit. 

Given the geometry of the sensor package, it is straightforward to calculate the portion of 
the celestial sphere covered by each sensor and the variety of independent measurements 
available for attitude determination. In particular, for the sensor package that is flown, 
there are two sun angle measurements over 50 percent of the sky and one sun angle measure- 
ment over the remaining 50 percent, such that the sun is fully covered with the package 
flown. In the slit sensor package, there are three sun angle measurements available over 87 
percent of the sky, one measurement over 12 percent, and there is no sun angle measure- 
ment at all for 1 percent of the sky. 

If the attitude is at orbit normal, and if the satellite is in an equatorial orbit, then in the 
package flown there are four sun-to-earth midscan dihedral angles available, and there are 
three such angles available with the slit sensor package. However, those three angles with 
the slit sensor package are available if it is required that the same sensor be used for both sun 
triggering and earth triggering. If you allow the sun to trigger with one sensor and the 
earth to trigger with another, then there are a total of nine sun-to-earth dihedral angles with 
the slit sensor package. 

There are two earth-width dihedral angles over 9 percent of the sky with the sensor package 
flown: There is one measure over just slightly less than half of the sky, and there is no 
earth-width measurement at all over 42 percent. On the other hand, with the slit sensor 
package, there are three earth-width measurements over 93 percent of the sky, two over 6 
percent, and one measurement over the remaining I percent, such that the sky is fully 
covered. 

In addition to this, there are three midscan-to-midscan nadir angle measurements over 93 
percent of the sky with the slit sensor package and one measurement over 6 percent. This 
measurement is not available at all from the attitude package flown. 

Thus, the results here indicate that the slit sensor package has substantially more independent 
attitude measurements over more of the sky than the attitude package flown, even though 
the package flown has more than twice the number of sensors. From an analytic point of 
view, the slit sensor package should give better attitude results from this information. 

Another characteristic of importance is sensor redundancy. That is, what information is left 
if the single most critical sensor for a particular measurement or a particular region of the sky 
is lost. With the sensor package flown, the loss of one critical sensor would be a moderately 
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serious problem, since we are guaranleed coverage of tlie earth over only 9 peicent of the 
celestial sphere. 1 he sun presents fewer problems, since there we are guaranteed coverage 
over 50 percent of the sky. With the slit sensor package, the loss of one critical sensor 
would be essentially no problem, since one of the sensors was intended to be redundant in any 
case. The earth and sun sensing are 100 percent covered; the sun angle measurement is at 
least 87 percent covered. 

The loss of two critical sensors would be a very serious problem for the attitude package 
flown. All of the sun or all of the earth observations could be lost. It would also be a 
problem for the slit sensor package. All of the sun angle and midscan-to midscan nadir angle 
measurements would be lost. However, the sun-to-earth dihedral angle and the e-- th-widih 
angles would still be available over at least 87 percent of the sky. Thus, it shocld still be 
possible to do attitude determination. However, the accuracy or bias determinatic.. char- 
acteristics could be impaired. 

If we lose three critical sensors in the slit sensor package, which has only three sensors, 
the spacecraft is obviously fully blind. In the attitude package flown, the spacecraft would 
be essentially blind. In some configurations, it is possible that there would be one observa- 
tion still in existence, however, it is unlikely that attitude determination could be done with 
that one observation. 

Lastly, there is the question of attitude accuracy, which requires a good knowledge of the 
sensor biases. A major factor in the accuracy with which biases can be determined is the 
coverage of the celestial sphere for each individual sensor. In the normal course of a mission, 
there are transfer orbits, inversion maneuvers, and so on. Each sensor will encounter a 
variety of geometries. As more coverage of the sky is available, more data from these 
different geometries will be provided, and the biases can be better determined, for two 
reasons. The most obvious, of course, is that there is more information available. In 
addition, there is a greater variety of geometries, and, in gen' ral, it is the variety of geo- 
metries which allows us to distinguish between sensor biases. 

Thus, it is worth examining the sky coverage of individual sensors for purposes of bias 
determination. With the attitude package that was flown on SMS, the sun sensors both 
covered 75 percent of the sky. When sensing the sun, the slit sensor package would cover 
anywhere from 87 to 100 percent of the sky. (The parallel sensor would cover 100 percent; 
the +8.5° sensor, 99 percent; and the -30° sensoi 87 percent.) Thus, the two packages 
would be essentially equivalent in this regard, with a very slight advantage, perhaps, going 
to the slit sensor package. 

However, when sensing the earth, there is a major difference. With the attitude package 
flown, the earth sensors cover only 20 to 29 percent of the sky. (The ±4° sensors would 
cover 29 percent; the +20° sensor, 28 percent; the -25° sensor, 27 percent; and the -48° sen- 
sor, 20 percent.) However, when sensing the earth, the slit sensor package would cover 93 
to 100 percent. With the parallel sensor, coverage would increase from 99 to 100 percent if 
an intensity measurement in the vicinity of the poles is used, rather than simply a pulse 
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measurement. So it is clear from this information that the slit sensor package would provide 
much better bias determination, for two reasons: There are fewer biases to be found, and 
those which are necessary are easier to determine. 

Thus, for SMS, it is seen that the slit sensor package has less than one-half the number of 
sensors that the flown package has, yet these would provide more independent attitude 
measurements over more of ihe sky, substantially greater redundancy, and far better bias 
determination characteristics than the flown sensors. One advantage of the slit sensor pack- 
age is that it is inherently simple. (It is possible to envision a prototype made out of a tuna- 
fish can and a photocell). A second major advantage is that it is exceptionally versatile in 
two senses: It obviously has full sky coverage, which allows for more independence of 
mission details, but probably more important is the great mechanical versatility. That is, 
liie same sensor can be used with changes in mechanical or electrical parameters to fit the 
particular mission at hand. Therefore, a single ground processing package could be used for 
a large number of missions. 

For example, we could change the tilt of the second sensor to provide the best pole-to-pole 
coverage for a particular mission. On some missions, it might be possible to use a tilted 
sensor that has two stops for different conditions where accurate attitudes are desired. Also, 
the triggering levels could be changed to fit the particular needs at hand. In theory, the 
same sensors aboard the same spacecraft could be used in transfer orbit 322 km (200 miles) 
above the surface of the .n as halfway to Mars, by simply changing the triggering levels 
and triggering on Mars and Jupiter as point sources and measuring their nadir angles with 
respect to the spin axis of the spacecraft. 

The slit sensor works on any spinning spacecraft. On any despun spacecraft, it can work 
by either rotating the sensors or by allowing them to look into one or more rotating mirrors. 

It is difficult to find any disadvantages from an analytic point of view, although, of course, 
such a sensor system may be physically unbuildable. One minor disadvantage is that a slit 
sensor’s response to the earth is not a square wave, as it is with point sensors. It rises to a 
maximum at the center of the earth and falls off toward the edges, which implies that the 
slit sensor might yield substantial triggenng level biases. However, triggering level biases are 
by far the easiest to resolve by ground processing, and therefore this is unlikely to be a 
serious difficulty. 

In summary, the optical slit sensor appears to be a versatile idea. It may be remarkably 
close to what is needed for modern space flight -a single standard attitude sensor with 
enormous versatility in its application. 
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ON THE DEVELOPMENT OF PRACTICAL NONLINEAR FILTERS 


H. W. Sorenson 
Applied Systems Corporation 
San Diego, California 


The general problem of estimating the state of a nonlinear, time-discrete system from noisy 
measurement data is considered from the point-of-v^w of developing feasible computational 
algorithms for evaluating the Bayesian recursion relations. Algorithms which have been 
proposed are reviewed, the computational implemtntaion of these algorithms is discussed, 
general conclusions coming from numerical studies are noted, and areas requiring additional 
research are defined. 
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CONVERGENCE CHARACTERi:>TICS OF BATCH AND 
SEQUENTIAL ESTIMATION ALGORITHMS 


B. Schutz 

University of Texas 
Austin. Texas 


The convergence rate and radius of convergence of the batch, and the extended sequential 
estimators are compared. The convened. :e behavior of the two processors in the presence 
of both observable and unobservab:*; parameters is considered. 
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MANEUVER STRATEGY DESIGN FOR MARINER/JUPITER/SATURN 
AUTONOMOUS GUIDANCE AND NAVIGATION 


T. Hagar 

Jet Propulsion Laboratories 
Ptsadena, California 


Candidate maneuver strategy algorithms for multiple and quasi-adaptive midcourse man- 
euvering are presented. Application of these techniques to the Mariner/Jupiter/Saturn 1977 
mission, and as possible candidates for autonomous navigation and guidance, is discussed. 
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CONSIDERATIONS FOR LARGE SPACE TELESCOPE 
(LST> MISSION EFFECTIVENFSS 

J. Tuttle 

' irlin Marietta Corporation 
Denver. Colorado 


In designing the support systems module for the LST, consideration must be given to hard- 
ware limitations and mission design requirements. Special software has been developed to 
analyze these sometimes conflicting requirements. Tliis talk will discuss the software and 
analysis made for the LST study. 
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SOLAR ELECTRIC PROPULSION 


Richard W. Barbieri 
Goddard Space Flight Center 
Greenbelt, Maryland 


Solar electric propulsion is certainly not a new concept. Indeed, it has been with us since the 
early 1900*s. But what is new is a growing awareness that solar electric propulsion offers a 
rather interesting alternative approach to study the earth and its environment and the solar 
system. 

This paper will cover some problems that we face in low-thrust mission analysis. After some 
preliminary comments about hardware and per jrmance parameters, concern will be devoted 
to the development of a nominal low-thrust trajectory and to the guidance and navigation 
problem. 

To put things into perspective, we should first discuss the major components of a solar 
electric propulsion system. It can be broken down into three major subsystems: One is a 
primary power source, which could be made up of batteries, solar cells, and reactors. Its 
function is to convert thermal or solar energy into electrical power. The second subsystem 
is a powe; conditioner and electrical control, which supplies specified levels of voltage and 
power to the heaters, valves, and thruster electrodes. In effect, it is an electrical power 
distribution center. The third subsystem is the engine, in which are mcluded the thruster, 
fuel tanks, and propellant and control system. 

Figure 1 shows a comparison of some performance parameters of both chemical and ion pro- 
pulsion systems. It can be seen that the chemical systems have a thrust-to-weight ratio rang- 
ing from about 10 to 10“^ g, operating with a specific i.i- ulse (Igp) in a range of perhaps 
90 to 250 seconds. The ion propulsion systems, on the other hand, operate with a thrust- 
to-weight ratio of about 10“^ to about 10*^ g, with a specific impulse ranging from 2000 
or 3000 seconds up to as high as about 1 1 ,000 seconds. 

The nuclear propulsion systems fit in this range of I to 10 g. Arc jets also fit into this range 
and overlap the chemical systems and ion propulsion systems as far as thrust-to-weight ratio 
is concerned. 

The low thrust system then provides a very high total velocity increment, and it does this by 
providing AV at a very low acceleration over a long period of time at very high specific 
impulse. The high specific impulse in effect translates to a smaller amount of propellant 
that has to be carried on board. 
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Figure 1 . Comparison of performance parameters 
of chemical and ion propulsion systems. 

Figure 2 shows what might be expected from a low-thrust propulsion system. The terminal 
mass-to-initial mass ratio is plotted against flight time in days; parameters are given for 
specific impulse and also for input power to the thruster-to-initial mass ratio (P/M^ ). 

The Solar Electric Rocket Test-C (SERT-C) mission is a study to place a spacecraft into a 
3!00-km circular orbit and slowly spiral out to geosynchronous orbit. The spacecraft will 
lift off with roughly 821 kg (1810 lb), with a specific impulse of roughly 3000 seconds 
for the transfer orbit engines, which arc about 30-millipound thrusters, and with an expected 
flight time of about 290 days. The SERT-C has a P/M^ ratio of roughly 4.2 and a termini,! 
mass-to-initial mass ratio of about 0.85. It will get into synchronous orbit with roughly 
703 kg (1550 lb) after lifting off with about 821 kg (1810 lb), a fairly high payload ratio. 

At the beginning of prelaunch analysis is the task of generating a nominal trajectory, which 
is usually optimal in som*' sense. This is where we encounter our first set of problems. 

There are certain phenomena peculiar to low-thrust problems, which must be modeled if we 
are to simulate a low-thrust trajectory with any semblance of accuracy: The first three 
items— geopotential, N-body, and solar radiation pressure— are not really peculiar to low- 
thrust systems but certainly must be included in a nominal trajectory algorithm. We have 
some experience with ballistic high-thrust-type missions with these three items. 

The next item, solar array interactions with the environment, are peculiar to low-thrust 
systems, since the array is the source of power to all spacecraft systems. A model of the 
radiation belt is required here, and the development of such a model will be akin to the 
development of the atmospheric density models we have had over the last 1 5 years. In 
both cases, we try to construct models for stochastic processes. Atmo.«pheric scientists may 
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Figure 2. Parameters of (ow-thrust propulsion system. 

have ex| .rtise in the development of radiation belt models, but those of us who are con- 
cerned with guidance and control and spacecraft systems are relatively unfamiliar with these 
models: t or instance, what kind of assumptions must be made to develop a working 
model to be inserted into a trajectory generator algorithm. 

Anotii . r phenomenon to be considered is the solar array degradiation, which is closely re- 
lated to the radiation belt model. For example, silicon solar cells can withstand a radiation 
dosage of perhaps 10'^ 1-MeV electrons. However, the radiation that might be expeiienced 
is perhaps two orders of magnitude larger than that. Such a dose, 10'* 1-MeV electrons, 
will have dire consequences on the available power to the thrust ~s. In particular, it could 
degrade the solar array power by as much as 40 to 50 percent of its beginning-of-life power, 
which, for the SERT-C mission, for example, is about 9 kW. Consequently, we are forced 
to protect the . olar cells with a thin coating of material about 76 to 152 nm (0.003 to 0.006 
inch) thick. Even so, the degradation must be modeled and will be a strong function of the 
thickness of the protective coating and of the type of material used in this coating. 

The last item to be considered is shadowing, which has been encountered before with regard 
to solar radiation pressure. No we must be concerned about it to determine the solar array 
and thruster performance during passage through shadow regions, in particular, during the 
thruster on/off times. The question is how long it takes the systen" to get up to full power 
after passage through shadow, and it happens that we do no‘ ^ave an answer at this time. 
Upon exit from shadow, for example, we know that it is going to take at least 10 mmutes 
to go thro’ -’h a preheat and controlled loop sequence. It is after this time that the thruster 
’•”ii operate at full power, which is expected to be a function of duration in shadow. 

. r -^f’cin here is the behavior of the a priori thruster biases after restarting the 
• :*' ’ ons are that such biases remain the same after restart, but detailed inves- 

' < M*'' warranted. The implication is that if such biases do change, then, for 
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the earth orbiting mission with solar occultation and thruster restarts quite frequent, the 
orbit determination process must reestimate these biases. Such frequent estimations could 
lead to significant orbit uncertainty, and this, in turn, has serious implications on the 
guidance policy. 

Having made these comments about the trajectory generation problem, we now turn our 
attention to the guidance and navigation aspects. I think it is safe to say that, of the small 
amount of work that has been done in the past in low-thrust mission analysis, little has been 
devoted to guidance and navigation. The problems here are difficult, and the opportunities 
for optimization studies abound. 

Environment and degradation models have already been mentioned. The same models which 
reside in the trajectory generator algorithm could certainly be used in the guidance and 
navigation algorithm. It must be emphasized that, in missions of this type, guidance and 
navigation are strongly coupled together because of the presence of a stochastic, continuously 
acting force. 

Another aspect of guidance and navigation is the thrusi vector model, which can be structured 
as a constant, as a constant plus noise, as a first- or second-order Markov process, or as a 
fully stochastic phenomenon. The first option is quite unrealistic. The fourth option leads 
us to extremely difficult mathematical problems, since it forces us to integrate random non- 
linear differential equations— nonlinear differential equations are difficult enough. Thrust 
magnitude depends on ion beam current, total accelerating potential, mass utilization 
efficiency, effective specific impulse, and numerous other parameters, each possessing a bias 
and time-varying components, which, when combined, may yield a standard deviation of 
about 5 percent of nominal thrust magnitude. But this is just a preliminary e^^timale. It 
could possibly get worse than 5 percent. The pointing error, on the other hand, a function 
of launch vibration, thermal distortion of the grids, and accelerated grid wear, in addition 
to improper knowledge or measurement of the thruster misalignments, gyro drift, misalign- 
ment, and other parameters. Therefore, modeling the thrust vector as a constant is not 
realistic. 

The next item is thruster orientation with respect to the solar array. If the engines are not 
gimballed, then the optimal orientation of the solar panels will not induce an optimum 
orientation of the thrust vector an^’ v^ce versa. This problem is one that really warrants 
many trade-off studies. Even if engines are gimballed and some freedom for the thrus- 
ters is allowed, optimization and trade-off studies must be carried out with respect to the 
relative orientation of the thrusters with respect to the solar array at various points of the 
mission. 

The last point to be discussed is the type and number of observations, two factors strongly 
affecting the navigation accuracy, which significantly impacts the guidance policy. One 
particular data type is that obtained from accelerometers. Because a low-thrust vehicle is 
thrusting over long periods of time, and because small deviations in the thrust direction and 
magnitude significantly alter the trajectory over these long periods, it becomes important 
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to evaluate the influence such data have on navigation error Implied here is another 
problem: When extremely sensitive accelerometers (sensing 10'*® to 10'** g with 2-arc- 
second accuracy) become (light-ready, a heavy burden is going to be placed upon attitude 
control sensor accuracy and measurement process. 

If only one accelerometer is placed on board along the nominal thrust axis, then information 
about mass flow rate becomes available (provided thrust magnitude is known fairly well): 
however, no information about thrust misalignment is available. On the other hand, it is 
expected that the navigation problems can be alleviated somewhat by placement of three 
highly sensitive accelerometers on board to reduce thrust vector misalignments. Sucii an 
alleviation is contingent upon precision alignment with respect to attitude control sensors 
Thrust direction with respect to the accelerometer axis can then be accurately determined, 
depending upon accelerometer accuracy, and referenced to inertial coordinates by the 
attitude control system. This is an area where very little work has been done and numerous 
studies must be made using not only earth-based data but also onboard navigation sensors. 

I would like to close by saying that the overall problem of low-thrust mission analysis is 
quite fascinating with new and nontrivial aspects requiring the development of new tech- 
nology. This is an area where it is necessaiy to reconsider a lot of the concepts that might 
have been formed in studying ballistic high-thrust-type mission analysis. The problems 
are not insurmountable, but they are going to be very time-consuming to overcome. 

This paper has discussed some of the mathematical models which are needed. In addition, 
there is the orbit determination problem where data types must be evaluated and used in 
combination with an optimal filter. Deciding upon a particular filter is not as easy as it 
might seem, taking into account the thrust-vector-related biases and time-varying compon- 
ents that must be estimated. The strong coupling between navigation and guidance and the 
pioblems it poses to the attitude control system arc crucial. 

DISCUSSION 

VOICE: What kind of funding is avail'»ble to study these types o( problems? They sound 
very interesting. 

BARBIERI: The funding right now is nebulous, at best. At this time we are not quite sure 
where we stand with regard to funding for low-thrust mission analysis. 

VOICE: Is anyone in particular interested in pushing this concept further? 

HOUGIIY: Yes. There is a possibility of a new start for solar electric research and develop- 
ment in the 1976 budget, but there is a very low probability of it actually coming into being. 
If it does not happen in the 1976 new start, it will probably be continued as a low level 
technology-type effort. The primary person who would be supporting it. should it stay at 
lew level technology, would be Jim Lazar at the Office of Aercnautics and Space Technology. 
From our point of view, there is a wide degree of uncertainty as to where wc arc going right 
now. We have to wait until we get a reading from the Office of Management and Budget 
on how they feel about it. 
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LOW THRUST OPTIMAL GUIDANCE FOR GEOCENTRIC MISSIONS 


T. Edeibaum and S. W. Sheppard 

Massachusetts Institute of Technology Charles Stark Drapper Laboratories 

Cambridge, Massachusetts 


Low thrust propulsion appears to have useful application as a means of satellite maneuvering 
in a strong gravity field. This thesis investigates the usefulness of one possible guidance 
scheme for such applications by means of a computer simulation. The guidance scheme uses 
some of the recent optimal trajectory theory applied to a particular class of orbit transfers. 
These transfers, between inclined circular orbits, are considered because they typify many 
mission objectives and have a relatively simple optimal solution. The optimal solution is 
presented here along with a mathematical approach to solving it on a computer The simula- 
tion program, which investigates the effects of an oblate gravity field on the guidance, is also 
presented. However, oblateness was found to cause relatively small errors and “closed-loop” 
guidance offered no significant improvement over “open-loop.” 
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RECENT INTERPLANETARY LOW THRUST STUDIES AT AMA 

F. I. Mann 

Analytical Mechanics Associates. Inc. 

Seabrook, Maryland 


Performance character!, tics of optimal low thrust rendezvous missions to the comets 
Giacobini-Zinner, Borrelly, and Tempel (2) with launches in the 1981-1986 time period 
are discussed. 

Also discussed are performance characteristics of optimal low thrust extra-ecliptic missions, 
including launch declination effects and the importance of optimizing the launch date. 
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GEOMETRIES DESCRIBING AN ORBITER’S RELATIVE MOTION 

J. B. Eades, Jr. 

Analytical Mechanics Associates, Inc. 

Seabrook, Maryland 


Analytical solutions to a set of moditied Euler-Hili equations lead to interesting geometric 
descriptions of a relative motion. Traces on the displacement and hodograph planes, de- 
fining a time history of the motion state, tell much of what can be expected from the solu- 
tion to any relative motion problem. 

The neoclassic solution of Clohessy and Wiltshire (for intercept) has been extended to include 
effects of forces and general initial values. These results are depicted on both the “local 
rotating” frame of reference and the companion “inertially oriented” one. 

General results for the relative motion state will be described, some special cases will be 
noted, and examples of uses of tf ese results will be mentioned. 
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STABILITY OF RELATIVE MOTION 


V. Szebehely 

University of Texas 
Austin, Texas 


The equations of the relative motion of two bodies in a given force-field are formulated, and 
it is shown that the conventional methods of representation lead to inst^ibility at rendezvous 
in the Earth's gravitational field. A method for selecting new dependent and independent 
variables is offered in order to stabilize the equations of relative motion. 
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LONG PERIOD NODAL MOTION OF SUN SYNCHRONOUS ORBITS 


Kenneth I. Duck 

Goddard Space Flight Center 
Greenbelt. Maryland 


The -jun synchronous orbit been used since the early 1960s for almost every meteorologi- 
cal satellite launched as well as for the two Earth Resources Technology Satellites 
(Landsat-I and Landsat-2). This well known orbit concept makes use of the earth's oblate- 
ness to induce a precession of the orbit line-of-nodes in order to maintain a fixed angular 
' .entation of the orbit plane relative to the mean sun. The sun synchronous orbit, even 
though it has been used for many years, will continue to be used for most earth observation 
and remote sensing applications because: 

• The sateUite passes through each latitude point at the same local time thus 
ensuring similar ground lighting conditions on each pass and consequently facili- 
tates data comparisons. 

• The average spacecraft solar array angle of incidence to the sun remains within a 
fixed boundary thus ensuring the availability of electric power. 

• The orbits (altitude and inclination) can be chosen such that the majority of the 
earth's surface can be mapped with near north-south contiguous swaths in a fixed 
period with repeatability. 

The long-period perturbations which disturb sun synchronous orbits have not been the 
subject of detailed investigations in comparison to geosynchronous orbits where the litera- 
ture abounds with material. The disturbances acting on sun synchronous orbits has not been 
examined in detail for several reasons: 

• With the exception of the Landsat satellite, previous sun synchronous vehicles 
have not had orbit adjustment capability to fight the perturbations. 

• The spacecraft hardware design lifetimes ( I to 2 years) have been sufficiently 
short that the perturbations could be neglected without significantly degrading 
mission performance. 

• Prior to the use of inertial guidance on the Delta launch vehicle, injection uncer- 
tainties were sufficiently large so as to mask any perturbation effects present. 

The assessment of the perturbations acting on sun synchronous orbits becomes more signifi- 
cant when longer lifetime spacecraft are developed as anticipated over the next 10 to 20 years. 


ca 
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With improvements in hardware technology and reliability and the cm-orbit refurbishment 
capability associated with the advent of Space Shuttle, operational satellite lifetimes could 
easily exceed 5 years. This increase in lifetime is almost a reality today as most earth 
orbiting spacecraft launched in recent years have exceeded their design life. 

The object of the study documented- here was to determine which perturbations significantly 
affected the long term nodal motion of sun synchronous orbits and then construct an approx- 
imate model which described the phenomena observed. Many computer simulations were 
made with several independent computer programs to assess the relative effect of various 
combinations of po'turbations. Typical of the perturbations included in the simulations were 
zonal and tesseral gravitational harmonics, third-body gravitational disturbances induced by 
the sun and moon, aixl atmospheric drag, it was observed that a model consisting of even- 
zonal harmonics through order 4 and solar gravity dominated the nodal motion. It was further 
observed that for long runs the orbit inclination and orientation of the line-of-nodes ex- 
hibited an oscillating behavior each having the same period. For all the cases run, the inclin- 
ation amplitude was very small (always less than 1 degree); however, the nodal motion could 
be quite large. Due to these observations, it was felt that a resonance existed between the 
inclination and the nodal motion. The mean daily rate of chaitge in inclination due to solar 
gravity (in radians/mean solar day) was found by analytic averaging to be 

di vl , 

— = 16200 — (1 cos i )* sin i sin 2 
dT 1} * “ 

where 

= mean motion of earth about the sun 

1 } = mean motion of satellite about the earth 

ij = solar obliquity 

= o’clock angle, the angle between the longitude of the ascending node and th: mean 
sun used computing local time 

i = orbit inclination 

The precession rate of the line-of-nodes due to zonal harmonics through order 4 (assuming 
a circular orbit) is 



where 



and 

Jj , = zonal gravitatkHuI coefficients 

Rc = mean equatorial radius of the earth 

a = orbit semimajor axis 

Differentiating $2, aibstituting di/dt, and noting that 
where 0.9856 “/day, then 

- - 

= - 16200 — (1 + cosi )* sin* i(A -3Bcos* i) sin2 
where A, B are A and B in radians/mean solar day. 

Examination of shows that for altitude regions where drag can be neglected, and 
assuming that sin* i is approximately constant, 

fioc h sin 2 

a form of the familiar pendulum equation. It is known that systems which are characterized 
by the above equation can exhibit libration or circulatory characteristics. Both characteris- 
tics can be observed on the phase plane plot (figure 1) made for ITOS-type orbits. It is 
observed that this system has stable equilibria which correspond to orbits whose line-of-nodes 
lie through the 6:00 a.m. and 6:00 p.m. points. The pendulum equation can be solved 
analytically using elliptic functions with a typical solution fca- a quarter cycle oscillation 
being 
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where 


m = 1/2 (1 + cos 2 12 ) 

o <T<P^/4 
o <n <90’ 

oc 

Note that P,^ denotes the libration period and denotes the value of o’clock angle where 
S2^ = 0. The libration period, is 

4K(m) 


where K (m) denotes a complete elliptic integral of the first kind. Figure 2 shows the 
o’clock angle libration period for ITOS-type orbits (h = 1489 km, i = lOI .9°). it is seen that 
the libration period increases from 26 years as the reference o’clock angle moves from one 
of the stable equilibria toward one of the unstable ones, investigation of the libration 
period has shown that the minimum libration period lies between 22 to 30 years for sun 
synchronous orbits between 200 and 2000 km. 


O 
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£ - 0.2 
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Figure 1 . O'clock angle motion in phase space for ITOS type orbit. 
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Figure 2. Libration period for ITOS-type orbit. 

The pendulum analogy has been compared with both simulation and flight spaceaaft data. 
Figure 3 compares the approximate solution with several simulations having various combin- 
ations of perturbations. The reference orbit in this comparison is the ITOS-type (h = 1489 km, 
i = 101 .9°) having an ascending nodal crossing at 3:00 p.m. local time. It is seen that there 
is excellent agreement between the pendulum analogy and those simulations which include 
zonal harmonics and solar gravity. 

The remaining figures (4 through 9) compare the o’clock angle time history generated using 
the pendulum model with that obtained from Brouwer mean elements for several flight 
spacecraft. In examining these comparisons one observes that there are two curves repre- 
senting approximate nodal drift propagation. The dashed curve is the propagation calculated 
using a set of orbit elements at the initial epoch for the particular satellite. The solid 
curve uses an iterated value of inclination to improve the agreement. This approach was 
taken because of the large relative uncertainty in measuring the inclination. It is seen in all 
cases that excellent agreement is obtained. In examining these data, note that the ESSA-2 
solution (figure 4) is a case where the motion lies in the circulatory region in phase space 
while the other solutions lie in the libration region. The ESSA-8 solution (figure 7) exhibits 
the turnaround expected for the libration motion. The final observation to be noted is for 
Nimbus-5 (figure 9) where the node is located near one of the unstable equilibria. There is 
still good agreement for this case where the disturbing acceleration is near zero. 

In conclusion, the nodal motion of sun synchronous orbits has been investigated and found 
to exhibit the characteristics of a pendulum. This pendulum motion results from solar 
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gravity inducing an inclination oscUlation which couples into the nodal precession induced 
by the earth's oblateness. The pendulum model has been compared with simulations and 
flight data with excellent correlation observed. 


AmtOXlMATC SOiUTKM 


A ESMATS»MULATION/ZOMALS.SlJN,MOGM 
O ESMAP SmULATION/ZCMAtS. SUN 
□ ESMAP SMULATION/ZONAU. MOON 



Figures. Approximate solution. 


EENOULUM ATPROXIMATiON 

(NOT DtSCERNIBLE ON THIS PLOT) 



Figure 4. ESSA-2 pendulum solution Gomparison. 
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Figure 5. Tiros-M pendulum solution comparison. 


KMOULUM AmtOXMlAriOM 

KMOULUMAmiOXtMATION 

WITH ITCRATCO IKlTIAL iWCUHATION 

A BASED ON ESSA-7 FLIQHT ELEMENTS 



Figure 6. ESSA7 p idulum solution comparison. 
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PENOULJM APPROXIMATION 



Figure 9. Nimbus-5 pendulum solution comparison. 


90 



N76-10176 


COMPARISON THEOREMS. NUMERICAL INTEGRATION, AND SATELLITE ORBITS 


Arnold Stokes 

Georgetown University 
Washington, D.C. 


Consider the equation 


x = X(x). (1) 

Assume that the behavior of the family of solutions is known, expressed by the estimate 

lx(t,x^)-x(t,x,) I < Kq f(t) lx, -Xg I , t >0, 

Here K„ is constant, and f is normalized by f(0) = 1 , f (t) > 1 . The function f (t) then 
expresses the stability or instability of ( 1 ), so that stability is equivalent to f = 1 , and a 
typical instability would be f (t) = t + I . 

Suppose we wish to numerically compute a particular solution of (1), say x (t, x^), for some 
fixed x„ . A remarkable theorem of Babuska asserts (under reasonable assumptions on X) 
that if a strongly stable difference scheme is used, numerical integration of (1) with x (0) = x„ 
is equivalent to obtaining an exact (theoretical) solution of 

y = X(y) + G(t,y), (21 

where y (0) = x, , (here lx, - x„ I is the starting error), and I G (t, y) I < rj, (here r) (small) 
reflects the discretization and truncation error). 

Babuska used this result to obtain global error bounds for (2), assuming x (t, x„) is a 
uniformly asymptotically stable soluvion of (1). Here we wish to compare solutions of (1) 
with solutions of (2), with particular attention given to the effect of different types of 
f(t). 

To develop the comparison, a converse theorem using Lyapunov functions due to Yoshizawa 
and Hale is applied. Let z = x - x (t, x^ ), then ( 1 ) becomes 

z = Z (t, z), (3) 

and the new estimate becomes 

• z (t, t^, z„) - z (t, t„, z,) I < K„ f (t) I z, - Zo' 
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f(0 , 

= K (tg) I z, - z^ I exp (a (t) - a (t^)) 


( 4 ) 


where K (t^) = f (t^), a (t) = log f (t). 

Then one can show there exists a function V (t, z) satisfying 

I zl < V (t, z) < K (t) I zl . 

I V (t, z) - V (t, y)l < K (t) I z - yl , and 


(a) 

(b) 


V 




(c) 


where Vj denotes the derivative of V along solutions of (3). Note that (a) and (c) allow 
one to recover (4), using standard theorems. Further, it is easy to show, using (b) that 


Vj <iv + K(t)n= Y 

Then an immediate consequence is that 

I y (t, X,) - X (t, x„)l < K„ f (0 1^1 X, - x„l + T?tj 

So we see that the term Kj, f (t) I Xj - I reflects the propagation of the starting error, 
and the term r? f (t) gives an estimate of the effects of the discretization and truncation 
error. 

Evidently, it is better to integrate a stable system, f = 1, than an unstable system, 
f (t) = t + 1 . A simple example is the unperturbed two-body problem, which is unstable, 
as the difference between two periodic solutions of different periods grows (at least locally) 
like t + 1. An example of a stable problem is the stabilized Kepler problem, as given by 
Baumgarte and Stiefel, where now all solutions have the same period. 

To extend the above to perturbed problems, consider 

i = X (x) + eX, (t. X) (5) 

Assume that in some region R, I X, (t, x) I • ' M, Here one supposes this region R contains 
the solution of interest and is large enough so that continued integration would not be of 
interest long before the integrated solution leaves R. 
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This assumption will allow the effect of f (t) on the error to be observed, for now we 
consider 


y = X (y) + eX, (t, y) + G (t, y). (6) 

and y (0) = X, , I G (t, y) K t? as before. 

The same V-function leads to the estimate 

I y (t, X,) - X (t, Xj,)l < K„ f (t) 1 ^ I X, - x„l + t (tj t eM)j , 

which is valid as long as both solutions remain in R. Ar^ain one can see the different error 
estimates obtained for f = 1 o*" f = t + 1 . Remark ; The foregoing gives upper bounds on 
the error. The strength of the conclusion then rests on the sharpness of the upper bounds, 
for one cannot conclude Error, < Error^ on the basis that a crude upper bound of the first 
is less than a crude upper bound of the second. 

However, by considering simple examples, one can see that at least in these cases, numerical 
integration of a stable system gives an error growth 0 (rjt), while for an unstable problem, 
the error grows like 0 (rjt^ ). So to that extent the foregoing estimates seem reasonable. 

Further, numerical experiments are being conducted at Goddard Space Flight Center using 
the stabilization techniques of Baumgarte and Stiefel, and the results should give further 
insight into the significance of the above results. 

In this regard, perhaps it should be observed that a similar approach to obtain lower bounds 
does not seem possible, as the difference cf two almost periodic functions will, in general, 
have an arbitrarily small lower bound. 

Note that the above . .imates can also be obtained in more complicated problems, with 
such means as time-dependent potentials. Details will appear in a forthcoming paper. 
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OmMAL EXPLICIT RUNGE-KUTTA METHODS 

D. Bettis and D. HaO 

University of Texas 
Austin. Texas 


Optimal explicit Runge-Kutta methods are developed for solving the initial value problem 
for systems of ordinary differential equations. These methods have an optimal estimate of 
the local truncation error term, thereby allowing the option of implementing a variable- 
step strategy. The coefficients of the Rui^e^Kutta method are selected so that the local 
truncation error is minimized and that the absolute stability region is maximized. 
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STABILIZATION BY MODIFICATION OF THE LAGRANGIAN* 


Joachim W. Baumgarte 

Mechanik-Zentrum, Technische Universitat Braunschweig 
Braunschweig, Federal Republic of Germany and Swiss Federal Institute of Technology 

Zurich, Switzerland 


ABSTRACT 

In OTder to reduce the error growth during a nummcal integration, a method of stabilization 
of the differential equations of the Keplerian motion is offered. It is characterized by the 
use of the eccentric anomaly as an independent variable in such a way that the time trans- 
formation is given by a generalized Lagrange formalism. The control terms in the equations 
of motion obtained by this modified Lagrangian give immediately a completely Lyapunov- 
stable set of differential equations. In contrast to other publications, here the equation of 
time integration is modified by a control term which leads to an integral which defined the 
time element for the perturbed Keplerian motion. 

INTRODUCTION 

It is well known that the classical differential equations of the Keplerian motion are 
unstable in the sense of Lyapunov. In general, Lyapunov-unstable differential equations 
develop more unavoidable numerical errors during a numerical integration than Lyapunov- 
stable equations do. We consider here the stabilization of the differential equations of 
Keplerian motion with the aim of improving the accuracy and efficiency during the 
numerical integration. We propose a stabilization method which is purely conservative in 
contrast to other methods (Baumgarte, 1972a). It is characteristic for all conservative 
methods, that they make the revolution time independent of the initial conditions 
(Baumgarte, 1974). 

GENERALIZED LAGRANGIAN 

In this method the stabilization goes hand in hand with the introduction of a new independent 
variable s instead of the time t. This procedure is called time transformation and s is 
called fictitious time. Furthermore, we will require that our stabilized equations of motion 
be developed from the Lagrangian formalism. But, here we have to use an appropriately 


* This paper was supported by the National Research Co*jndi and the National Aeronautics and Space Administration. 
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modified generalized Lagrangian formalism. In order to introduce, instead of t, the 
fictitious time s as a new independent variable, we have to use, instead of the original 
Lagrangian L 


L = L (q., Qj, t). 


I = 1 , 2 ,. . . ,n. 


( 1 ) 


where q. are the coordinates and dot means differentiation with respect to t, the following 
generalized L^rangian L*: 


L* = ht' + p 




( 2 ) 


In equation (2) the prime means differentiation with respect to the independent variables. 
The time t is now a dependent variable (time coordinate). Together with t appears its 
conjugated momentum h, which represents physically the negative total energy. 

Furthermore, we consider only the case, where 

/i = /i(qj,h)>0 (3) 


is a freely chosen positive function only dependent on q. and h. Later we will see that /r 
may be interpreted as the local scale of the time transformation. More general dependences 
of the scale p are also of interest but shall not be considered further here. Our special 
choice of p retains the equivalence between Lagrangian and Hamiltonian formalism. Th' 
choice has two consequences: 

1. The transformed kinetic energy in an inertial system is also a quadratic form in 
the velocity components q' with respect to the fictitious time s. 

2. \ conservative system remains conservative after the time transformation. 


With these restrictions we obtain the differential equations of motion in the following fOTm: 


d 

ds 



dL* 


d 



dL* 


= 0 


= 0 


£ 

ds 





(4a) 

(4b) 

(4c) 


Relative to the system (4a, b, c) the following remarks are in order: 

• Equation (4a) states the equations of motion, which are also given by the 
original Lagrangian L (q., q^, t) after changing from t to s, but which are 
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distinguished only by control terms (Baumgarte, 1972b). This fact will be observed 
li'^er in the special case of the Kepler problem. 

© ' quation (4b) gives: 

h' = p(q.,h) ^ LCq-.ql.h,!), (5) 

dt 

which asserts the energy relation. In the conservative 
ase, where the original Lagrangian L does not depend explicitly 
•n the time t, we get h = constant. 

• From equation (4c) there follows: 

. ( dp d ) 


It will be obvious that the expression in the bracket is a control term, which 
represents the energy relation. This control term vanishes in the exact solution 
of the equations of motion, in such a way that t' = fi. Therefore, it follows as 
previously promised that is the local scale of the time transformation. 


EXAMPLE: KEPLERIAN MOTION 

We consider now the keplerian problem and use Cartesian coordinates. We will choose the 
scale in such a way that the fictitious time s will be for a pure Keplerian motion the 
eccentric anomaly. Only by doing this can we obtain both Lyapunov-stable differential 
equations and equivalence between Lagrange and Hamilton. We, therefore, choose: 


p=— ,r=lxl, (7) 

y/lh 

where r is the ciistanro 

Before we establt.. the generalized Lagrangian L*, we first write the original Lagrangian L. 
With the masa .n = I , K as the gravitational parameter, L has the form: 

I 

L = -lxP+— . (8) 


With the help of equation (2) we now obtain as the generalized Lagrangian L*: 

h 


, y/lh , ,2 

L* = ht + I X I + 

2r " 




(9) 
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With (9) the stabilized equations of motion are: 


(x.5) 



"] 


(lOa) 


h' =0 


(10b) 


5 




Ik^ 

It*' 



(lOc) 


(10a) requires that h* = 0. It is essential that in contrast to an earlier publication (Baumgarte, 
1972b) the vectOT-eq nation for x" does not depend on h. This means the revolution time 
is the fixed number 2rr and this implies Lyapunov-stability for the x" - equation. 

From ( 1 Ob) there follows h = constant. The constant h is computed from the initial con- 
ditions and is placed in the computer as a fixed value. This is the presupposition that now 
in contrast to earlier publications (Baumgarte, 1972b) the time integration (10c) is 
Lyapunov-stable. The expression: 

j_ ( .!£!!• 

2v^l ^ i 

a part of the right hand side of (10c), is proportional to the difference between potential 
and kinetic energy. Therefore, the mean value with reject to the fictitious time s of this 
difference vanishes because we have an oscillator problem in principle. Consequently, the 
expression K} /(2h)^^ in (10c) represents the exact mean value of the right hand side of 
(lOc). Tills fact implies Lyapunov-stability also for the time integration. 


(X)NTROL TERMS AND LYAPUNOV-STABILITY 


In order to show the effect of the stabilization by the control terms, in the classical 
Keplerian equations 


K* 


X = - 




i. 



( 11 ) 


we substitute in place of the independent variable t the fictitious time s by using 
dt/ds = t' = r/y/2h! Because h = constant we obtain: 


(x.x') 


x' + 


2h r 


0 


(12a) 


98 



h' = 0 


(12b) 


=0. (12c) 

V2h 

We now transpose system (10) in such a way that we can see clearly that system (12) differs 
from system ( 1 0) only by control terms which are the right hand sides of equations (13a) 
and (13c). 


(x-x*) 

X 

Ix'l^ 


1 


X + — 

J.2 ~ 2h 

r 

2r^ 

2hr 

+ — X 
2 ‘ 

(13a) 


h' = 0 




(13b) 

t'- — = 

Ix'l^ 


1 

r 



2r* 

" 2hr 

+ — “ 
2 


(13c) 


In the control terms in ( 1 3a) and (13c) the same bracket appears as a factor, which is 
analytically zero with respect to the energy relation. 

Equation (13c) or (10c) can be integrated. We find as an integral of motion: 



The control terms in (13a, c) produce the Lyapunov-stability under the supposition that 
the constant h is computed once and for all from the initial conditions. The proof for the 
stability of the complete system ( 1 3) or ( 1 0), respectively, (with respect to h = constant, 
whereby this constant is to be computed, finally, by the initial conditions) can easily be 
carried out by making the transformation into action and angle variables, because the 
equivalence between Lagrange and Hamilton exists. By doing this, the generalized 
Hamiltonian, obtained from the generalized Lagrangian by making a Legendre transforma- 
tion, will be linear in the action variables, which implies Lyapunov-stability (Baumgarte, 
1972b). Another proof follows from the equivalence of the system (10a, b) together with 
( 1 4), to the corresponding equations of the KS-transformation (Stiefel and Scheifele, 1971), 

We will call attention to the fact that the dependence of the time transformation 

t' = xfy/^ on the momentum h makes possible the elimination of the instability under the 

presupposition that h' = 0 is integrated exactly. 

In the case of perturbed Keplerian motion, the stabilized system ( 1 3) or ( 1 0) is modified 
by additional perturbation terms. In equation (1 4), C is no longer constant but becomes 
the slowly varying time element (Baumgarte, 1972b; Stiefel and Scheifele, 1971). 

99 



Numerical experiments have always shown a reduction in error in the numerical integration. 
It appears that the positive effects of the stabilization of the pure Kcplerian motion carry 
over to the perturbed problem. 

We will finally remark that the KS-transformation can be performed directly in the general- 
ized Lagrangiap L* by inserting 




(15) 


thereby giving immediately the analogous complete, stabilized, set of differential equations 
which leads directly to the KS-elements(Stiefel and Scheifele, 1971). 
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SUMMARY 

A set of initial Cartesian coordinates, which are free of ambiguities and resonance singularities, 
is developed to study satellite mission requirements and dispersions over lon£ lifetimes. 

The method outlined herein possesses two distinct advantages over most other averaging pro- 
cedures. First, the averaging is carried out numerically using Gaussian quadratures, thus 
avoiding tedious expansions and the resulting resonances for critical inclinations, etc. 

Secondly, by using the initial rectangular Cartesian coordinates, conventional, existing 
acceleration perturbation routines can be absorbed into the program without further mod- 
ifications, thus making the method easily adaptable to the addition of new perturbation 
effects. 

The averaged nonlinear differential equations are integrated by means of a Runge Kutta 
method. A typical step size of several orbits permits rapid integration of long lifetime 
orbits in a short computing time. 

INTRODUCTION 

Several sets of averaged elements (Lorell, 1970; Broucke and Cefola, 1972; Uphoff, 1973) 
are in use for satellite lifetime studies. These usually suffer from ambiguities and resonance 
singularities for low inclinations, near circular orbit, near polar orbits, critical inclination 
resonances, and such. Moreover, it is necessary to develop the perturbation representations 
in these element coordinate systems, which often requires ingenuity and is difficult in ap- 
plication. The Cartesian coordinates have been extensively utilized and routines are avail- 
able to numerically generate most of the significant perturbations. The initial conditions 
of the Cartesian solution of the classical two-body problem have been developed for variation 
of parameters (Pines, 1961 ; Christensen, 1970;Godal and Johansen, 1968), but has not 
found wide application. With the advent of averaging as a tool for eliminating long tedious 
numerical integrations in computing solutions, this study was undertaken to reestablish the 
initial Cartesian coordinates as a useful set of parameters for orbital analysis. 
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THE INITIAL CONDITION CARTESIAN ELEMENTS 

The equations of motion of the satellite in the planetary reference frame are given by 

R = -#i — + F (1) 

where F represents the perturbation forces other than the central attraction of the principal 
body. 

The initial Cartesian coordinate parameters which describe the motion are given in terms of 
the position and velocity vectors in Cartesian coordinates by 


Rq = -fR + fR 


( 2 ) 


where f, g, g, f are given as functions of the difference in eccentric anomaly, 0, as 


f =1- 


a(l - cos 9 ) 


1 ^ 

g =— (rv/^a sin 0 - d a(l - cos 0)) 
P 

V/Ja sin 0 


'o' 


g =1" 


a(l -COS0) 


d =R*R 


('■0 


d cos 0 


- a(l - COS 0) + r cos 0 - sin 0 

0 VM 


1 ' ^0 _ 2 R R 

a /i t ti 


(3) 


The relationship of the time increment to 0, in the absence of perturbations, is given by a 
form of Kepler’s equation 
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(4) 


t-t 


0 




r d 

0 -sind +— sin^ 

a 


(1 - cosfl) 


The differential equations for the variation of and under the action of F are generated 
in Pines (1961) from the conditions that 

d 

— R = 0 
dT 


and 



(5) 


Thus, we have 

^R„=g^R-g^R-gF 


d . . 

— R„=-fR + fR + fF 
dr ° ^ 


( 6 ) 


Following Godel and Johansen (1968), we choose for the perturbation equation of the 
difference in eccentric anomaly, 6 


6 


= 0 


(7) 


This serves to eliminate mixed secular terms (see Christensen, 1970) from the perturbation 
derivatives of f, g, f, g, that arise in Pines (1961) where the time from the initial time was 
assumed unperturbed. 




(8) 


We can replace the vectors R, R from (6) using the inverse of (2) 

R = fRo+gRo 

R = fRo+gRp 


(9) 


^ Rq = (g/ - fg,) Ro + (ggr ' 8r8) Ro ‘ 


^Ro=(f,^ff)Ro^(M-f^)Ro + fF 


(9a) 
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An important stability condition on (9a) is derived. We have 

fg-gf= 1 


It follows from the perturbation derivative of (10) that 

f^-f^ = -g/+fg. 


( 10 ) 


( 1 !) 


Thus the coefficient of in the d/dr expression is the negative of the coefficient of 
Rg in the d/dr R^, expression in Equations (9a). This will reduce the computing work in 
the averaging process. 


Since the perturbation of the transit time is not zero, we require a differential equation for 
the time. Using (3), (4), and (5), we obtain the perturbation differential equations for the 
time, 



1 

V^T 



- sin 0) + sin 0 - 

2>/r 


d(l -cos0)\ 

v#r 


a(l -cosd) 

Vir 


+ 1 


(12) 


where 


a. =2a^ R F 


= R • F 


The perturbation derivatives of f, g, f, and g are given below: 


% 


f^ = — (1 -COS0)I — 


^0 » 


/rsinO d(l-cosO)\ a(l -cos0) 
g=f )a_- — d 


(i 


r T 

IJi 


ra ‘0 
f =fj-L - — 
" Ua r. 


(1 -COS0) 

c = - a 

^ r 


(12a) 


(13) 
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where 


-cos0)a^- 


V^sinfl 

nAT 



(13a) 


We now proceed to the averaged differential equations. 

THE AVERAGED EQUATIONS OF MOTION 

The principle here is to replace the nonlinear differential equations for R^fT) and R^,(t) by 
their average value over a single period of the reference orbit as defined by Rq(t, ) and 

Let the period of i 'erence orbit be T(t, ), then for 


T T 


(14) 


Thus, 


i. “ 

dr 


t-t=r-T= — 


VH 


(6 - sin 0) + r sin 0 + 


d^a(l - cos^) 




j; / / T'"'’ 


d \ f” i 

dr^*"‘»'"27T j 


d0 


105 



where 


r *io 

— *I-cosO + — cos 9 -I- — sm( 
* * \'ia 


(15. 


Equation (15) couk) be integrated analyticaUy, using Fourier series in 9. However, this 
would require a representation of F in d and greatly burden the introduction of additional 
perturbations. Moreover, mathematical resonances would appearand require special 
techniques for each resonant term. In this, wc follow the lead of Uphoff ( 1973) and adapt 
Gaussian quadratures s the technique few evaluating the integrals. 

We note that F is a vector in three space, and recalling ( 1 1 ), we need to evaluate only 1 0 
int^rals. The average equations become: 

=a,Ro + a,R^+C 


d - 

— R =a,R„-a.R„+F 
dr 


3 ‘ 0 


(16) 


where 


-(.-to)=l.a. 


y*ir 

‘• ’h / 


/•* 

’^ 2 ,/ 




r. / 


(1 6a) 


G = 


-r. / 
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I 

2it 


1 


2n 




3 


U6a) 

(cont.) 


To obtain the averaged F, G vectors, it is necessary to specify the perturbations. For the 
purpose of this study, we consider accelerations due to tesseral and zonal harmonics for a 
rotating planet, third body forces, such as the sun and moon, and drag. The detailed 
equations are given in Appendix A. It should be noted that where F is a function of R, 

R and t, for the purposes of averaging, these functions are given by the two body equations 
and Kepler's equation given in (17a) through (I7d) referenced to the r time at the midpoint. 

We now consider the numerical integration of the averaged equations. We propose to use 
a Rui^e Kutta method with a r step equal to several periods. It must be borne in mind that 
the time variable is r and not time. For each evaluation of the seven derivatives ( 1 6), 

at a specific r time, Tj, we must carry out the averaging procedure. Using Rq(t.) R^(r.) as 
the reference orbit, we compute the 10 integrals given by (16a). 

The Runge Kutta solution will produce Rj(r), Rj,(r) and (t - 1^,) (t). To compute an 
ephemeris of the state, R (t) and R (t), we proceed as follows: 

Let the period of the R„(v), R^fr) orbit be T (t). 


where 


T(t)= — (a(T))3/^ 
P 



R„(r)R„(r) 


Let N be the integer part of t/T (t), and d is given by, 

= 2v N + a 

Using Newton’s method, solve Kepler’s equation for a 


At = t-NT = 


v5T 


(o - sin o) + r^a'^ sin o + 


d^a( I - cos a) 

y/ir 


(I7a) 


(17b) 


(17c) 
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For any position in the N + I revolution, where the incremental eccentric anomaly, 6, lies 
between 0 and 2 tt. 

Let 


t = t(T| + 


r(T) d(T)(l - cos a) 

o-sina+— sino+ (l+a^(r)) 


y/ia 


R(t)=f(«)R„(r) + g(a) R„(r) 


R(t) = f(a) R„(r) + g(tt)R„(r) 


(I7d) 
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APPENDIX A 


We consider several typical |x:rturbations. 

Tesseral and Zonal Harmonics. 

Let the central body rotate about its polar axis, k (t) with a uniform angular rate cj. If the 
inertial Cartesian coordinates of the vehicle is R, then the unit vector to the vehicle in planet 
fixed coordinates is given by, 

R 

u =k(t) — 
r 


R 

s =i(t)* — 
r 


(A.l) 


R 

w=j(0 — 
r 

where i (t) and j (t) are orthogonal reference axes fixed in the rotating body equation, 
perpendicular to the polar axis k (t), expressed in the inertial Cartesian system of the 
vehicle. 

Following Pines, ( 1973). the acceleration in the inertial system is given by 

R(t) 

F| =a,i(t)+ajj(t)+ajk(t)+a^ — — (A.2) 

For the purpose of completeness, we list the a. coefficients for J, . J , . . and . 

a,( ) = ttjO j ) “ ^ “ “il-* J ) ~ *1 1-*4 I = ^ ® 

ftal 3Jj (I - 5 U-) 




2r" 


a4(J3) = - 


fia^ 5Jj(Jtu-7u^) 


(A.3) 


a.f J = - 

4 4 


fal 15 (-1 + I4u^ -21 u'*) 


Sr'* 




pa^. 3JjU 
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(A3) 

(toat.) 


3Jj (5u^ - I) 
“3<J3> = -— ^ 


nil 5J^ (7 u* - 3 u( 

“3(^4) = -— ; 

2r* 

fiaj6 

&,(t Sjj) = — (wSjj + s^22^ 


paJ-6 


«,(C,,.S^) = 


(sS 


22 


■wCj,) 


a,(C^,.S^) = 0 


lial 15 

a^(Cjj.S^) = -— — -(C^(s^ -w^) + 2swSj2) 

Third Body Acceleration. 

Let the gravitational mass coefficient of third body be u.. Let the ephemeris of the third 
body, R (t), be given with respect to the central planet, about which the vehicle is orbiting, 
then the indirect acceleration is given by 



(A.4) 


where 


r . - K • • 

Cl Cl 


r = IR-R,;! 

VI Cl 


(A.5) 


Atmospheric Drag. 

Let p be the air density, given as a function of the vehicle position. R, and the time. t. The 
drag acceleration is given by 


F 


3 


-jP(R,t) 


vCnS . 

R 

m 


(A.6) 


no 



where 


V = IRI 

Cp ' coefficient 

\w = vehicle mass 
s = effective vehicle arag area. 
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NONLINEAR COUNTER EXAMPLE FOR BATCH 
AND EXTENDED SEQUENTIAL ESTIMATION ALGORITHMS 


B. T. Fang 

Wolf Research and Development Corporation 
Riva-dale, Maryland 


A simple nonlinear example is presented which shows the welldcnown iterated batch least- 
squares and extended sequential estimation algorithms may converge to different estimates. 
For this example one may even say the extended sequential algorithm converges to the 
“wrong” value. 
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ORBIT DETERMINATION IN THE PRESENCE OF 
UNCERTAIN ATMOSPHERIC DRAG 


B. Tapley, D. Dowd, and B. Schutz 
University of Texas 
Austin, Texas 


Uncertainties in the knowledge of the atmospheric density in the associated drag parameter 
constitute one of the primary limitations on the accuracy on which the orbits of near earth 
satellites can be determined and predicted. In most orbit determination programs, the effect 
of uncertainties in atmospheric drag are determined by adopting a standard atmosphere and 
estimating the drag parameters, 0 C^ A -j- m. 

However, for most missions, Cp and A vary and the standard atmosphere will contain errors. 
Each of these factors will lead to errors in the orbit determination and prediction operation. 
In this presentation, an approach for estimating the drag parameter, the effective satellite 
cross sectional area in the atmospheric density simultaneously with the satellite state, is 
described. 
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APPUCATIONS OF SATELLITE-TO-SATELUTE TRACKING 
TO ORBIT DETERMINATION AND GEOFOTENTIAL RECOVERY 

P. Aigentiero, R. Garza-RoUes, and M. O’DeD 
Goddard Space Flight Center 
Greenbelt, Maryland 


Recent simulations have demonstrated the applicability of satellite-to-satelUte tracking data 
to the related problems of orbit determination and geopotential recovery. Specifically, 
satellite-to-satellite tracking between an earth orbiting satellite and a satellite at geosyn- 
chronous altitude (36000 km) produces long continuous data arcs which are not available 
by means of ground-based tracking. This facility, in conjunction with correct estimation 
techniques, can yield exceptional orbit determination accuracy. The data type also has 
considerable applicability to geopotential determination when the low satellite has a high 
inclination. 

ORBIT DETERMINATION APPLICATIONS OF SATELLITE-TO-SATELUTE TRACKING 

Attention is focused on the difficult problem of determining GEOS-C altitude with an 
average accuracy of 1 m. Error sources considered were ground -station survey error, data 
bias, epoch state errors for the high and low satellite, errors in spherical harmonic coefficients 
of the geopotential field to degree, and order 8. Standard covariance analysis software was 
utilized to determine that survey error and data bias were insignificant error sources, but 
that uncertainty in relay satellite state caused radial errors of the order of 100 m. Geopo- 
tential uncertainty caused an average radial error of 6 m. When the high and the low satellites 
are estimated and other error sources ignored, the average radial error is 6.2 meters. To 
identify the geopotential terms to be estimated in order to satisfy the constraint of 1 meter 
altitude resolution, a recursive procedure is implemented. If N-dominant geopotential terms 
along with GEOS-C and ATS-6 satellite states are estimated, and if the 1 meter radial 
error constraint is not satisfied, the geopotential term from among tne unadjusted set which 
causes maximum radial aliasing is added to the estimated set of parameters. If the constraint 
is still unsatisfied, the recursive procedure continues. This recursive procedure has been 
automated within the covariance analysis software. The result of the procedure when applied 
to the GEOS-C orbit determination problem is that after 31 dominant geopotential terms 
are recursively identified and added to the estimated state along with GEOS-C and ATS-6 
state, an average radial error of 1 meter is achieved. 
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GEOPOTENTIAL RECOVERY APPLICATIONS 


During the GEOS-C mission it is planned to track GEOS-C from ATS-6 from widely separated 
geosynchronous positions of 94°W and 34°E. The resultant data set should be almost globally 
distributed. Simulations demonstrate that it is possible to estimate from this data set coeffi- 
cients of the spherical harmonic representation of the geopotential field to degree and order 
8 with an accuracy on an order of magnitude superior to that presently obtainable. 
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LARGE ANGLE SATELLITE ATTITUDE MANEUVERS 


John E. Cochran*' and John L. Junkinsf 
The University of Virginia 
Charlottesville, Virginia 


ABSTRAC: 

Two methods are proposed for performing large angle reorientation maneuvers. The first 
method is based upon Euler’s Rotation Theorem; an arbitrary reorientation is ideally accom- 
plished by rotating the .pacecraft about a line (the “Euler Axis" or “Principal Line”) which 
is fixed in both the body and in r-.*ace. This scheme has been found to be best suited for 
the case in which the initial and desired attitude states have small angular velocities. A de- 
tailed evaluation of the associated feedback control laws and sensitivity to disturbances has 
been carried out, assuming the control system to consist of four single-gimbal control moment 
gyros (CMGs). These results indicate that the proposed scheme is feasible with realistic 
physical constraints on the CMG torque source. The sev«.md scheme is more general in that 
a general class of transition trajectories is introduced which, in principle, allows transfer 
between arbitrary orientation and angular velocity states. The method generates transition 
maneuvers in which the uncontrolled (free) initial and final states are matched in orientation 
and angular velocity. The forced transition trajectory is obtained by using a weighted 
average of the unforced forward integration of the initial state and the unforced backward 
integration of the desired state. 

Our current effort is centered around practical validation of this second class of maneuvers. 

Of particular concern is enforcement of given control system constraints and methods for 
suboptimization by proper selection of maneuver initiation and termination times. Analo- 
gous reorientation strategies which force smooth transition in angular momentum and/or 
rotational energy are also under consideration. 

DISCUSSION 

Many spacecraft must perform one or more reorientations or attitude changes during their 
lifetimes. The ways in which these maneuvers are performed are obviously important from 
the standpoint of conservation of energy. However, often the optimality of a maneuver in 


*Vi$iting Associate Professor, Associate Professor, Aerospace Engineering Department, Auburn University, Auburn. Alabama, 
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terms of the energy required to perform it may be of lesser importance than the time and 
computational power needed to define it. That is, often, finding the optimal maneuver 
strategy may not be desirable or even possible within the constraints imposed. 

In this paper, two methods for defining “good,” nominal, lai^ge-angle attitude aaneuver^ 
for spacecraft are presented. Nei"' er of the methods is generally optimal (although in 
special cases they would be), but both offer the advantages of being ( 1 ) relatively easy to 
use and (2) explicit, rather than iterative. The first method, which is well-suited for the 
case of quiescent initial and final rotational states, is based on Euler’s Rotation Theorem, 
that is, it is a single axis maneuver. The second method, which may be used when the 
initial and/or final rotational states ai'e not quiescent, is based on the us^ of transition 
trajectories in a phase space of dimension eight, wheie the redundant dimensions are due 
to the choice of four parameter descriptions of orientation. 

SINGLE AXIS MANEUVERS 

The idea of utilizing single axis rotations for arbitrary reorientations is not a new idea (Meyer, 
1966). For example, the Apollo Command and Service Module were reoriented with a 
single rotation about the required Euler axis (Crisp et ah, 1967). Such m.' neuvers are not 
necessarily optimal, in fact, Dixon et ai. (1970) have shown that single axis maneuvers of 
axisymmetric spacecraft through the use of thrusters are generally more costly in terms of 
fuel used than two-impulse maneuvers designed to minimize fuel expended. They are, 
on the other hand, more easy to define than optimal maneuvers, and if the spacecraft is 
asymmetric, no closed-form optimal maneuver strategy comparable to that of Dix* . et nl. 
(1970) is available. 

One important concept, used to some extent in both methods, is that of Euler, or eigenaxis, 
rotations. Figure 1 is included to illustrate this concept. On the left hand side of figure 1 
are shown a centroidal body-fixed system Cxyz and the associated unit vector triad (u, , u. , 
Uj) as well as the Euler axis for a partiailar reorientation. The Cxyz system pictured on the 
right-hand side is a rotationally inertial coordinate system which is used along with the 
displaced Cxyz system to indicate how the Euler axis rotation proceeds. 

Figure 2 provides some information concerning how the attitude of a moving trihedral 
C Cj e^ ej with respect to a fixed trihedral C e' e^ e^ can be defined. The rotation matrix 
A, whose elements are direction cosines relating the two trihedrals, can be constructed using 
Euler angles, say 6, and <p or alternatively using Euler parameters, here indicated by 
ttj, , «! , , and ttj . Furthermore, the Euler parameters are intimately related to the 

principal rotation angle 0 and the direction cosines fi, , Ej , and Cj of the principal line, i.e., 
the Euler axis (or equivalently, the eigenaxis corresponding to the unit eigenvalue of ^), 

Assuming that A is specified, the four parameters a^,, , and can be determined in 

a noniterative fashion. Then, for single axis rotations, the direction of the Euler axis can be 
determined as well as the required principal angle. For control of single axis rotations in 
the presence of disturbing torques, it is advantageous to define an eigenaxis system as 
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Euler Axis ^ 

or the 
Elqenaxis 

Figure 1. The Euler axis arrd Euler axis rotatioirs. 
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EigenaxTs 
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Figure 2 . Attitude change logic. 



deso'ibed in Cochran et al., 1 975. The way in which such a system is defined is briefly 
summarized in figure 3. 



Figure 3. The eigenaxis system. 

Figure 4 indicates how attitude errors in the form of small Euler angles and may be 
specified by using the body-fixed eigenaxis system Ce^e,^ e^ which ideally has its e^-axis 
directed along the principal line. 


tar.$„ - 2(a^a, + - aj * \) 

tan«^ * ' “m^ 

6 , 4 - "sman" error angles 
m n 



Figure 4. Attitude errors. 

The dynamical equations for a spacecraft which contains n control moment gyros (CMGs) 
are given in figure S. This equation, the control torque equation, and CMC steering law 
also given in figure 5 were used in a recent study (Cochran ct al., 1975) of the feasibility 
of controlling single axis rotations using CMGs and several linear controllers. 
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Figure 5. Attitude dynamics and CMG steering law. 

The following equations are the nonlinear matrix forms of the equations of rotational motion 
and approximate linear equations derived assuming small angular velocity magnitude, small 
errors and a linear contioller of angular acceleration (for more details see Cochran et al., 197S). 

s 0 

" hlf - !« Ic 

o' = 

’£ = I->T® + I-^T® 

^ =e -ex = -c 



Linear Feedback 



External Disturbance 

One of the controllers used in the study reported in Cochran et al. (1975) is depicted 
schematically in figure 6. The controller is a proportional-plus-integral-plus-derivative, or 
PID, controller modified by using a “model follower” commanded rotation angle generator 
which serves two purposes. First, it allows maintaining the difference in the actual angle of 
rotation 0 about the principal line close to the commanded rotation angle 0j (hence higher 
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gains and tighter control) and second, it provides a means of specifying the nominal rotation 
rate oi* a priori so that allowable CMG gimbal rates will hopefully not be exceeded. Atten- 
tion should be drawn to the form of the function . This function is such that the com- 
manded rotation angle connects the initial value for i.e., zero, and the desired final value 
of , with a smooth curve. 



Fixed PIO Corifiouratlon 

'•gure 6. PID model-follower. 


Some typical response curves * x the rotation angle are shown in figure 7. Four different 
controllers were used to generate the curves by numerically simulating the attitude motions 
of a Large Space Telescope type spacecraft using four single gimbal CMGs a.s torque sources. 
The full nonlinear equations were used in the simulations (the reader is again referred to 
Cochran et al. (1975) for more details). The PID model-following controller has been dis- 
cussed and the PID and PD are simpler controllers. In figure 7, MRAS refers to a model 
reference adaptive system controller which utilizes Lyapunov’s second method to generate 
variable control gains. 

In the study reported by Cochran et al. ( 1 975), the control of single axis rotation in the 
presence of disturbances was found to be feasible using rate limited CMGs. 

TRANSITION TRAJECTORIES 

The basic idea for the second method of performing attitude maneuvers was motivated by 
the use of the function in the PID model-following controller and previous use of the 
concept of weighting functions in the areas of geodesy and gravity modeling (Jancaitis and 
Junkins, 1974). In figure 8, this idea, the use of an averaging concept for definition of 
transition maneuvers, is summarized. The functions a, (t) and (t), shown in the upper 
left-hand portion of figure 8, represent an attitude variable, a (t), say one of the j = 0, 

1 , 2, 3, as it would appear as a function of time if the spacecraft were rotating freely in its 
initial state (subscript 0 and similarly its final, or desired, state (subscript b). The dashed 
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curve represents a transition curve for the variable a. Note that a (t) is constructed by 
using the unforced forward integration of a from time t^ using the actual value of o at t^ 
as an initial condition aitd the unfOTced backward integration of a from time t^ with the 
desired final value of a as a final condition. 




Figure?. Rotational response. 



The weiqht functions ouerantee satisfaction of 
the 2nd order osculation constraints: 


at t = t. 

at t = 

n(t«) » 

o(t,) • ^(t,) 

4!| 

dt|t, dt |t. 

!!2| - 

dt[t^ dt |t^ 


<l*3[ ^ '*'^1 

dt^t, dt* |t, 

dt'itf dt' It, 

for arbitrary qiven functions a^(t) and 


a^(t) * Unforced forward fnteqration of actual initial state at time t« . 
a^(t) • Unforced backward inteqration of desired flrwil state at time t^. 

..(t) “ "Transition Trajectory” from a^{t) to 


Figure 8. Averaging concept for definition of transition maneuvers. 
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The weighting functions, and , shown as functions of normalized time r, are used in 
defining a (t), so that the second-order osculation constraints (see right-hand side of figure 
8) of agreement in value and the first two time derivatives at times t^^ and t^ are satisfied. 

Note that W^. and are relatively simple functions and that they provide a method for 
defining a smooth function which generates a transition curve between two arbitrary, but 
smooth, functions of time. 

Recalling that a is an attitude variable (it may also be a vector of variables), it is obvious that 
the proposed method provides a means for explicitly determining the torque required to 
make the spacecraft attitude follow a prescribed trajectory in phase space. Thus, the trajec- 
tory departs smoothly from the initial state and ends in the desired final state. The explicit 
determination of the torque depends only on the availability of the first two time derivatives 
of three, or more, attitude variables. 

We assume that, during the time required for the maneuvers, the external tm-ques acting on 
the spacecraft are neglible. Furthermore, in obtaining the results presented here, we have 
assumed that the spacecraft is a single asymmetric rigid body with known ine.'tia character- 
istics. These assumptions allov. the use of an analytical solution for the free rotational 
motion of an asymmetric rigid body (see Morton et al., 1973, and Kraige and Junkins, 1974) 
in computing the forward and backward states and the necessary time derivatives used in 
constructing the torque required to perform the maneuver. 

The analytical solution for the torque-free rotational motion of asymmetric ri^d body which 
was used is summarized in the following list. 

Analytical Solution for Torque-Free Motion of Asymmetric Body. 


Anquiar Velocity 


Euler Orientation Parameters 


u * flt + 

Wi(t) = Wi,«dn(u,k) 
to»2(t) » W2RiSn(u,k) 
wi(t) - 


{a(t)} H 


where 



Yti 

Yi 

Y2 

Ys 


“Yl -Y2 -yT 



Y» “Y> Ya 



Y» Yo -Yl 



-Ya Yl Yft 




where 

are constants defined In terms of 
inertias and initial conditions. 


6o(t) ^ [1 ♦ asn(u,k)] * cos«o(t) 

Bi(t) 

6a(t) 

B,(t) 


^ [1 - asn(u.k)] cosai(t) 


% 1 

[1 ♦ asn(u,k)] sinto(t) 


j [1 - asn(u.k) j sin^, (t) 


ta ,Yo ,Yi »Ya »Y ^ } are constants defined by 
initial conditions and inertias. 


and 

♦o(t) and ♦i(t) are functions involving in- 
complete elliptic integrals of the third kind. 
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The body-fixed components of angular velocity are denoted by , and cOj and these 

are functions of the time t and constants defined in terms of principal moments of inertia and 
initial conditions. The Jacobian elliptic functions dn (u, k), sn (u, k), and cn (u, k) are the 
basic time functions involved in the solutions for u>, , , and cj, . 

For the attitude description, Euler orientation parameters were chosen. The body’s inertial 
attitude is defined by the oc^, j = 0, 1 , 2, 3, which are expressible as functions of a set of four 
constant Euler parameters 7j,j = 0, 1, 2, 3, (which define the orientation of a nonrotating 
angular momentum coordinate frame), and four time varying Euler parameters j = 0, 1, 

2, 3, (which define the orientation of the principal axes of the body relative to the angular 
momentum frame). The 0 . may be expressed as functions of time by using Jacobian elliptic 
functions, incomplete elliptic integrals of the third kind and, of course, initial conditions and 
inertias. 

Transition trajectories may, in principle, be defined in terms of any set of attitude variables. 
Two particular four parameter sets were chosen. The Euler parameters ot, j = 0, 1,2,3, 
are a rather obvious choice; however, preliminary studies have indicated that the set composed 
of the principal angle of rotation 0 and the directions cosines Sj, j = 1 , 1 , 2, 3, of the Euler 
axis may be a more desirable set. These will be referred to as the principal rotation paiam- 
eters. The construction of transition trajectmes using these two sets of variables is sum- 
marized in the following list. 

Definition of Transition Trajectories 

Principal Rotation Coordinates 
(ftot''i,ai2,a3) (iff 1 s) Transfomation 
J^(t) = 2cos*‘{a^^) I 

^ ^ 

(t) = /sir, ^ I 
Forward Trajectory 

*f{t) = Function (inertias, actual initial state at to. t) 

(t) = Function (inertias, actual initial state at to, t) 

' i 

Backward Trajectory 

t|j(t) = Function (inertias, desired final state at t^, t) 

(t) = Function (inertias, desired final state at t^, t) 

Transition Trajectory 

t(t) = W^(T)i^(t) + W^(T)t|j(t) 

i = 2. 3 

= w^(x)f^ (t) + W^(T)tjj (t) 

e.(t) = t/T- [t^(t) + 


i = 1. 2. 3 


i = 1. 2. 3 
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Euler Parameters 


Forward Trajectory 

(oi^(t)} = Function (inertias* actual initial state at to. t) 
Backward Trajectory 

{oijj(t)} = Function (Inertias* desired final state at t^, t) 
Transition Trajectory 



«o(t) = ±/l - [aj(t) + a,^(t) + a3(tn 


Note that when four parameters are used appropriate constraints must be introduced. 

The transition trajectory concept has been used to generate transition trajectories and the 
associated torque time histories for several pairs of initial and final states. Figures 9 and 10 
show results for a case in which the principal rotation parameters were used to generate 
the transition trajectory. The principal moments of inertia picked for the example are I, = 3, 
Ij = 2 , and I 3 = I • The maneuver required was to change the rotational state of the space- 
craft from one in which =v/3/2,«j = -y/l5/8, a = l/3,a^ = 0, CO, = 1, co^ = 0.10 and 
cOj = 0.0 at t = t„ = 0 to one in which = 1/4, <*, =s/lSI4,a^ =a^ =0,o3^ =0.101, 
cOj = 0.0 and cOj = 0.0 at t = t^ = 5.0. In this maneuver the inertial components of angular 
momentum were changed from (2.93, 0.47, - 0.48) at t = 0 to (0.303, 0.0, 0.0) at t = 5. 
Basically, a state in which the spacecraft is rotating rapidly with an initial nutation angle of 
about 1 2 ° is changed by the maneuver into a state of much slower (an order of magnitude 
less) spin about the principal axis of maximum moment of inertia, with the orientation of 
the bcxly’s angular momentum vector also changed. 

The time histories of the angular momentum components and the Euler parameters are 
presented in figure 9. The principal rotation parameters were averaged to produce the 
transition maneuver, but the angular velocity and Euler parameters are, of course, also 
averaged in the sense that the desired final state is reached in a smooth manner. On each 
plot in figure 9, the solid curve represents the time history of the indicated variable which 
results from unforced forward integration, the dashed line has a similar interpretation, but 
is derived from the backward integration of the desired final state, and the bold curve is the 
transition curve. Very smooth transitions of all seven variables are evident. 
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Figure 9. A - Angular velocity; B - Euler parameters. 


Time histories for the variables which were averaged (i.e., the principal rotation parameters) 
are shown in flgure 10. Also in figure 10 are included the time histories of the body-fixed 
components of the torque needed to generate the tranation maneuver. The three curves 
in each of the plots indicated as LHATl, LHAT2, LHAT3, and PHI are analogous ti those 
previously described. In the torque plot, the 1 -component is indicated by the solid line, the 
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2-componcnt by the dashed line, and the 3-component by the bold line. The difference in 
the magnitudes of the initial and final angular velocities is apparent from the radically 
different slopes of the curves for the forwaid and backward solutions for the principal angle. 



(A) 

Figure 10. A - Euler axis direction cosines; B - principal angle; C - torque history. 


Considering the torque history in more detail, we note that the second-order osculation con- 
straints embodied in the functions and result in zero torque at t = 0 and t = 5. All 
three components of the torque are smooth and bounded ; the largest magnitude of any one 
component is about 2. 1 5. 


127 



Referring to the list of transition trajectories, the transformation from the Euler parameters 
to principal rotation parameters embodies a singularity whenever 4> is an integer multiple 
of 2 jt; resulting in ambiguous definition of the Euler axis. Thus, in the absence of some 
remedial action, the transformation equations and their derivatives do not provide an accept- 
able basis for defining multirevolution transition maneuvers. We are currently studying means 
of circumventing this difficulty and thereby allowing application of this method to the 
multirevolution case. 

From the work which we have done to date on the attitude maneuver problem, we have 
drawn the following conclusions: 

• Single-axis maneuvers are not necessarily optimal, but provide a reasonable basis 

for quiescent-state-to-quiescent-state attitude maneuvers using onboard computations 
and continuous torques, especially if the spacecraft is asymmetric. 

• Control of such single-axis maneuvers in the presence of disturbances is feasible. 

• Transition maneuvo'S provide an explicit solution to a more general class of man- 
euver problems. 

• Control of transition maneuvers looks feasible. 

• An important feature of both methods is that iterative solution of a two-point 
boundary value problem (TPBVP) is avoided. 

• Transition maneuvers provide starting solutions (which satisfy the boundary 
conditions) for iterative solution of TPBVPs. 
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MISSION OPERATIONS FOR THE LOW COST 
MODULAR SPACECRAFT 


R. L. des Jardins 
Goddard Space Flight Center 
Green belt, Maryland 


The low cost modular spacecraft (LCMS) was developed by Goddard Space Flight Center to 
provide a standard spacecraft bus which could easily be conflgured to support virtually any 
near-space unmanned mission to be flown in the 1980s. The I.CMS features subsystem 
modularity allowing great flexibility and on-orbit servicing, yet achieving benefits of wide- 
spread standardization. Also, the LCMS design incorporates a high-performance onboard 
computer as a remote controller for most spacecraft subsystems. The LCMS is described 
briefly and its implications for mission operations is explored. 
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ONBOARD ORBIT DETERMINATION USING SERIES SOLUTIONS 


T. Feagin 

University of Tennessee Space lu -titute 
Tullahoma, Tennessee 


An iterative, linear sequential estimation algorithm is presented which is suitable for use in 
a small onboard digital computer. The solution is obtained in the form of a finite series of 
Chebyshev polynomials. The series-solution provides a close approximation to the actual 
orbit which is vaiid for a given interval of time. A Kalman filter is used to combine new 
observational data with the old estimate of the state and its associated error covariance 
matrix in order to update the series-solution and thus to provide a new optimal estimate 
and covariance matrix. 



HIGH ALTITUDE AUTONOMOUS NAVIGATION 


Howard A. Garda 
Martin Marietta Corporation 
Denver, Colorado 


PURPOSE OF STUOV 

The applied research described in this paper pertains to a high altitude autonomous na 'na- 
tion system which was the subject of a Phase 0 Preliminary Design and Feasibility Study 
under contract with SAMSO, United States Air Force Systems Command. This contract 
had the expressed objectives of selecting a particular design configuration by means of a trade 
study involving several candidate sextant concepts and carrying out a preliminary design based 
on the Anal selected sextant subsystem. In addition, the contract also called for an analytic 
evaluation of the general navigation concept by a numerical analysis which included a 
parametric sensitivity study and a performance demonstration by a Monte Carlo analysis. 

SYSTEM REQUIREMENTS 

Air Force requirements speciAed that the system be autonomous to the extent that no 
dedicated earth emission would be necessarily invoked and that the system would operate 
effectively at least 1 20,000 n.m. The space sextant system that was adopted relies upon 
no earth refer ;nced missions, either passive or dedicated. The navigation accuracy has been 
shown to be constant for any of the tested orbits, irrespective of orbit shape, orientation, 
or altitude. The Air Force requirements further speciAed that the system should have a 
demonstrated accuracy (by numerical analysis) of at least I n.m. ( I o rss) at a confidence 
level of 95 percent and, where large trajectory errors are presumed to exist, should converge 
in at least 10 hours from the onset of navigation. The numerical analysis has shown that 
the system converges to less than I n.m. in about 6 minutes and has demonstrated steady- 
state navigation accuracies of about 800 feet ( I a rss) after I S to 18 hours of measurement 
processing. Other design goals established for the system include mission versatility, 
satellite versatility, insensitivity to reasonable parameter variations, insensitivity to 
satellite maneuvers, a 5-year lifetime, and the utilization of existing technologies. 

NAVIGATION CONCEFI 

The navigation concept is simple and direct. Navigation is accomplished by means of 
me<. ared angles between the brighter stars (visual magnitude < 2) and the bright limb of 
the moon. Reduction to the moon’s center, including compensation for asphericity effects 
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and lunar terrain, is accomplished by onboard software. The essential data required to 
determine the navigation position are the measured angle, the moon’s ephemeris, and 
precise time. In principle, it may be shown that angular measurements from each of two 
stars to its nearest limb on the moon establishes a line of position for the spacecraft. 

Similar measurements made on the earth’s limb provide a second line of position. The 
intersection of the two obtains a complete navigation fix in as short a period of time as is 
necessary to complete these measurements. 

High navigation accuracy is achieved by further improving this position (and velocity) 
knowledge by recursively filtering subsequent star-moon measurements over the next 
several hours. 

ANGLE MEASUREMENT SUBSYSTEM 

The bask xtant instrument consists of two Cassegrainian telescopes, an angle measurement 
head, and two gimbals providing for two additional degrees of angular freedom. 

The electronics package consists of an oscillator, registers, A to D converters, a digital 
microprocessor, and the wheel speed control servo. The total device, sextant and electron- 
ics, would weigh less than 25 pounds. 

The principle of operation is quite simple. A spinning element inside of the measurement 
head, running at constant angular velocity, intercepts the optical path of a ray originating 
at a star observed by one telescope, and subsequently intercepts the path of a ray originat- 
ing on the moon’s limb observed by the second telescope. The time elapsed between the 
reception of these two signals, which may be recorded with great precision, is directly 
translatable into arc measure. iTiis associated servo system incorporates two independently 
operated subsystems: an in-plane servo which positions the star (S) and the limb (L) 
trackers precisely on these respective targets and a cross-plane servo which orients the 
measurement head (wheel and both optical tracking telescopes) into the piano of meas- 
urement defined by the S X L vector. 

MEASURING HEAD OPTICAL SCHEMATIC 

Two optical trains are utilized in the operation of each tracker, the tracker ray and the 
timing pulse ray. The tracking ray enters the telescope aperture, reflects off the primary 
mirror, then the secondary mirror, and finally impinges on the detector. The timing pulse 
ray, originating at an inteinai light source, passes through the collimator lens and is reflected 
off two mirrored surfaces in a prism that is common to both trackers and that rotates with 
the wheel. This timing pulse ray is then reflected off this primary and secondary and 
impinges on the same detector as the tracking ray. The detectors arc two-stage, four 
quadrant, differential detection types, with the first stage for acquisition and the second 
stage for ptecise tracking the timing pulse generation. The timing pulse is generated by 
zero crossiiifc detection of the timing pulse ray as it crosses two detectors whose output 
is differenced. One advantage of this type of detection is the accuracy that can be 
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achieved with zero crossing detectors and another is its insensitivity to a mismatch in 
detector output responsivity. 

The sextant trackers are designed for total symmetry and reciprocity so that either tracker 
can be used for star or limb tracking. This feature also allows for attitude measurement 
and onboard self'calibration using two stars. 

SYSTEMATIC ERROR COMPENSATION 

High measurement precision and stability is achieved by means of a phase-locked loop and 
self-compensation for radial runout and encoder disk systematic errors. The phase locked 
loop drives the wheel assembly at an angular velocity of approximately SO rad/s. The 
commanded rate originates with the oscillator and the position feedback is derived from an 
optical transducer disk by means of a read head. When the wheel is in motion, the output from 
the read head is frequency. The actual disk contains inscribed sin 2^^ and cos 2^^ functions. 
These functions permit the extraction of phase information and make a wideband, high 
gain phase loop possible. 

Systematic errors in the encoder disk are on the order of 10 arc-seconds until self-compen- 
sation procedures are activated. Other error sources which include bearing noise, sensor 
noise, and data sampling become dominant after self-compensation; however, the combined 
effect is to yield a total measurement error of about 0.5 arc second over a one second 
measurement time interval. The automatic tracking provision allows for continuous meas- 
urements during low thrust maneuvers. During high acceleration maneuvers, however, the 
sextant will be caged, but may resume operation within one minute following the maneuver. 

NAVIGATION PERFORMANCE 

The principal outcome of the analytic investigation of the system performance was to 
demonstrate that the system is capable of exceeding Air Force requirements by a sub- 
stantial margin, both with regard to navigation accuracy as well as convergence time. A 
Monte Carlo simulation consisting of 59 samples (with no outliers) was necessary to ensure 
a 95 percent confidence level in the order statistics, constituting the primary investigative 
tool used in the performance analysis. 

Two orbits were used as the models for the Monte Carlo analysis. Orbit A is a highly 
eccentric, 12-hour Molnyia type orbit, and Orbit B is an equatorial, circular, geosynchronous 
0 . 1 ^. t. The Monte Carlo analysis considered only star/lunar limb measurements; consequently, 
these results show the high dependence of convergence time upon orbit geometry, for instance, 
7 hours for Orbit A and 1 2 hours for Orbit B. By definition, convergence time is the time 
from the onset of navigation to the last major inflection of the navigation error curve; 
however, system accuracy continues to improve with additional data processing so that 
ultimate steady-state accuracy is finally achieved after several days of sextant navigation 
where the moon has completed a significant segment of its orbit. The accuracy is on the 
order of 800 feet and is a reflection primarily of the uncertainty in the lunar ephemeris 
itself. 
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A later analysis which employed star/earth limb measurements to augment the star/lunar 
limb measurements showed that the convergence time could be reduced to a few minutes. 
The final system accuracy would be achieved in less than 20 hours, owing to the vastly 
improved geometry resulting from the intersection of the second line of position. Ir 
either case, the ultimate navigation accuracy depends upon the mwe precise lunar Umb 
data which is the [wimary mode of navigation, the earth limb data being used only to speed 
up convergence and to supplement the lunar limb data when the moon is occulted by 
either the earth or the sun. 

STATISTICAL AGREEMENT 

The performance analysis used sev^aJ statistical methodologies to represent the expected 
system accuracy, all of which demonstrated good relative agreement, indicative of a 
fundamental internal consistency and providing a firm basis for the main conclusioi's drawn 
from the study. Four curves, ail representative of the I a case, were computed and plotted 
for the two demonstration orbits: 

• 67% - 67% population coverage derived from histogram data across 

59 samples; 

• o - ensemble statistics standard deviation across 59 samples; 

• - square root of the summed position eigenvalues from a covariance 

analysis; and 

• O 3 - square root of the second central moment derived by a time averaged 

(moving window) technique along a single sample simulation. 

An apparent disagreement between the covariance analysis and the Monte Carlo results at 
the initial time was due to the methods used to initialize the respective trajectory errors. 

The initial state covariance for these orbits was set at the estimated 1 a values ; whereas, the 
initial trajectory offset for simulation purposes was set at 3 o values for all six state elements 
(a probability of occurrence o.* only 7 * Kf*®). Furthermore, the curve, being one 
sample from the Monte Carlo set, is also shown to deviate initially from the other curves. 

The important result, however, is the fact that all of the statistical methods converge to 
essentially the same values at steady state for each of the two performance test cases. 

Orbits A and B. 

SENSITIVITY ANALYSIS 

A sensitivity analyas was carried out on three high altitude orbits (Orbit C, 14,350 n.m.; 
Orbit D, 68,000 n.m.; Orbit E, 1 15,000 n.m.) in addition to Orbits A and B. First and 
foremost, these sensitivity analyses have demonstrated that the system is totally indifferent 
to altitude in the earth-moon domain. Of course, this result is not unexpected because 
the prime observable is the moon and not landmarks on the earth. In the absence of aug- 
mented observations using the earth’s limb, the rate of convergence appears to obey a 
logarithmic function based on orbit period. However by exploiting earth limb measurements. 
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atl orbits converge in the period of a few minutes and obtain steady-state navigation 
accuracies of less than 0.2 n.m. in a matter of hours. 

The senytivity analysis also considered a range of model parameter errors of up to twice 
nominal values. These parameters included the lunar ephemeris ( I o = 600 feet), lunar 
terrain height, ( 1 o = 1 700 feet), sensor noise ( I o = 0.566 arc sec), and initial trajectory 
errors, the latter being dependent upon the particular orbit type. The principal results 
were as follows: 

• The system is virtually insensitive to expected lunar ephemeris errors, even for 
the twice nominal case; 

• 2 X r .iinar terrain errors are acceptable, but this also means that onboard 
compensation for terrain height will be required; 

• Nominal sensor errors arc near optimal in view of the contribution of other 
modeled system errors; and 

• Large dispersions in initialization errors were completely suppressed by the 
recursive filtering process. 

OCCULTATION OF THE MOON BY THE EARTH OR SUN 

A phenomenon that must be acknowledged by the system concerns the possibility of an 
occultation of the moon by either the earth or the sun. The effect upon the navigation 
system was investigated by simulation, and it was found that the navigation accuracy was 
largely unaffected so long as the augmented measurement mode (earth limb) was employed 
during these critical periods. 

In the specific example tested, worst case geometry was assumed for Orbit B (24-hour, 
circular orbit), where the earth itself is twice occulted by the sun in 24 hours, and the moon 
is concurrently obscured by the sun for 30 hours. Prior to either of the occultations in 
this example, the large initial trajectory errors were reduced by combined earth limb and 
lunar limb measurements. The earth and moon were then presumed to enter simultaneously 
the 9° look angle constraint zone centered about the sun. No navigation data could be 
acquired until the earth was again visible after 1.5 hours, and a slight increase in navigation 
error may be noted. Earth measurements were then commenced for the next 22.5 hours 
when it was again occulted; however, the growth of navigation errors this time were 
sufficiently suppressed so that no noticeable increase in navigation error occurs during the 
second occultation. Lunar limb measurements were resumed after 30 hours. 

PRECISION ATTITUDE REFERENCE SYSTEM 

The principal feature of the Space Sextant Navigation System is the high precision of the 
angular measurement. This same feature conduces to make the device an attitude 
reference sensor as well. The latter capability shows promise of becoming one of the most 
accurate onboard attitude sensing systems in existence. As an attitude sensor, it will be 
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necessary to add a platform reference mirror and a precision three axes gyro package. An 
autocoiiimation light source would also be added to the basic instrument in order to permit 
one telescope at a time to align itself with the reference mirror. Two axis angular reference 
may be achieved either with a second orthogonal reference mirror or a precision base encoder 
to measure the yaw angle. The relative advantages of these two alternative modes of 
obtaining two angular measurements of a single star have not been assessed at this time; 
however, it is anticipated that the system will ultimately be capable of 0.1 arc sec, three 
axes orientation under steady-state conditions. 

SS-HANS DESIGN SUMMARY AND CONCLUSIONS 

A fundamental included angle accuracy of 0.5 arc sec (1 a) by the Space Sextant makes 
navigation practical in cislunar space to a high degree at accuracy (< 0.2 n.m.) and at the 
same time provide the highest level of autonomy possible. An attitude reference system 
of 0.1 arc sec precision is also feasible, employing the basic sextant instrument modified 
to perform autocoiiimation in conjunction with a fixed mirror system. The sextant and its 
associated electronics will weigh less than 25 pounds and wUl require 7.5 watts average 
power, 50 watts peak power, and 30 watts during actual measurement. Design lifetime is 
5 years with redundancy and parts derating, bearings included. 

The Space Sextant is designed to be functionally subservient to the spacecraft computer 
system and to impose no cycle time restrictions on the computer. Software functions 
including navigation, on-orbit calibration and health monitoring will require 15K of 16 bit 
words of memory. Read-only mass storage for a 5-year lunar ephemeris and lunar terrain 
height data will require up to 140K of 30 bit words. The latter needs may be fulfilled 
either by a solid-state ROM of less than 4 pounds or by magnetic tape units. The read-only 
provision greatly enhances the long life and reliability of either type of mass storage device. 
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UTILIZATION OF LANDMARK DATA IN 
ATTITUDE/ORBIT DETERMINATION 


N76-X0181 


H. Siddafingaiah 
P. S. Desai 

Computer Sciences Corporation 
Silver Spring, Mary bind 


Picture data are taken both by geostationary and medium altitude satellites primarily for 
meteorological purposes. Abundant availability of picture data suggests that, both for nav- 
igation and picture registration for earth resources and meteorological purposes, the 
possibility of using landmarks is attractive. (Techniques of pattern recognition for identify- 
ing known landmarks are being developed currently at Goddard and by other contracting 
companies and universities, therefore, my talk will not address itself to this aspect.) One 
may find usage of landmarks in: 

• Reduction of tracking requirements, 

• Autonomous navigation, 

• Potential improvemeni in attitude accuracy, and 

• Improvement in geographic analysis of picture data. 

One has the following options in using picture data from satellite fixed camera to determine: 

1 . Position given attitude, 

2. Attitude given position, 

3. Both attitude and position given: 

(a) enough landmarks 

0?) realistic attitude and orbit dynamic models 

The problem is to express landmark coordinates in terms of satellite position, attitude, and 
camera scan angle. 

This presentation will describe a mathematical model for determination of satellite position, 
velocity, and attitude using landmark coordinates as observables. This model, although 
developed with respect to earth-stabilized missions, Tiros-N and Nimbus-G in particular, 
is applicable to any earth-stabilized satellite in general. 
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As is usual in any mathematical discussion of flight mechanics, one has to use several coordin- * 
ate systems simultaneously. The coordinate systems used here are: 

1. (x‘,y',z‘) - Inertial, earth-centered 

2. (x'^ , y*^ , z** ) - Rotating, earth-centered 

3. (x® , y® , z® ) - Satdlite body fixed 

4. tx^ , y^ , z^ ) - Instantaneous auxiliary orbital, earth-centered 

5. (x®,y®,z®) - Instantaneous orbital, earth-centered, 

where 



If n is a unit vector along the line of sight of the landmark, then, 

n® =G(o,0,7)n® (1) 

where G (a, P, y) is the attitude matrix associated with 2-1-3 rotation. 

If T), and f are the angles made by the line of sight with the positive direction of the axes, 
then 



( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


(7) 
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we have 
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where 


Finally 
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where 

E is the 3-1-3 rotation matrix 

R is the transformation matrix from the inertial to the earth-fixed system 

0 is the longitude of the ascending node, and 
ir 

^ = w + — 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 
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where 


cj is the angle of perigee 

V is the true anomaly 

In the light of relations (3) and (6), one can also consider the camera angles as observables. 
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ABSTRACT 

Given the locations of several landmarks on a satellite acquired image and their true geo- 
graphic coordinates, the position and orientation of the satellite can be determined. Two 
methods for automatically locating the image coordinates of specified landmarks are des- 
cribed. The first, a particularly fast sequential similarity detection algorithm for template 
matching was originally described by Nagel and Rosenfeld. The second method involves 
iteratively resampling the picture function in the vicinity of the anticipated landmark. A 
variety of other speedup methods is also described. An application to SMS imagery is 
envisioned. 

INTRODUCTION 

An important role in satellite navigation is played by the recognition and location of land- 
marks on images of the earth’s surface as viewed by the satellite. Given the locations of 
several landmarks on the image, and their true geographical locations, the position and 
orientation of the satellite can be determined. The mathematics of this process, and its 
error analysis, will not be recapitulated here. (See, for example, Phillips and Smith, 1972.) 

We shall assume here that the approximate position and orientation of the satellite are 
already known, and that we want to use the landmark data to obtain more accurate estimates. 
This implies that we know approximately where, on the image, landmarks should appear- 
perhaps to within a few dozen picture elements (pixels). We shall ignore here any errors 
in our estimation of the orientations and sizes of the features; it is reasonable to assume 
that if the error in estimated position is small, then the errors in estimated orientation and 
size are negligible. 

If we regard the orientations and sizes of the landmarks on the image as shown, and their 
positions as approximately known, then the problem of landmark identification and loca- 
tion can be solved by a template matching process. This basically involves comparing a 
template or small image of the landmark with the picture, in the range of positions where 
the landmark could be located. If a good match is obtained, the landmark has been detected, 
and the position of this good match gives the location of the landmark. 
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There are many different template matching processes that could be used for this purpose, 
and the matching can be implemented in various ways (digitally, optically, or by a human 
operator). We shall consider here only automated digital techniques. 

The objective of this proposed effort is to make the template matching as computationally 
efficient as possible, while keeping machine storage requirements at the minicompuier level. 

The standard digital template matching technique, which is based on cross-correlation, is 
relatively slow. Substantially faster results can be obtained by using an absolute difference 
measure of mismatch (sometimes known as the sequential similarity detection algorithm 
(SSDA). This approach is not only cheaper computationally, but also lends itself to the 
u^ of efficient search techniques for speeding up the detection of mismatches, as shown 
by Nagel and Rosenfeld, 1972. It can be combined with various pre- and post-processing 
techniques that should increase the sharpness of the matches obtained, thus increasing the 
accuracy, as well as the speed, with which landmarks can be located. 

TEMPLATE MATCHING 


Use of the Absolute Difference Misnnatch Measure 

The most commonly used measure of the match between a template, t (x, y), and a oicture 
p (x, y), is the correlation coefficient (or normalized cross-correction) 



(x, y) p (x + u, y + v) dxdy 



(x, y) dxdy 



(x + u, y + v) dxdy 


where (u, v) is the displacement of the template relative to the picture. Here the area of 
integration, T, is the area occupied by the template. This coefficient takes on values be- 
tween 0 and 1 , and has value 1 only whe^^ t (x, y) and p (x + u, y + v) are identical (over 
the area T) except for a multiplicative constant. In the digital case, the expression for the 
correlation coefficient is identical, except that the integrals are replaced by sums over pic 
ture element arrays. 

Computation of the corr :!at*'on coefficient is relatively slow, because it involves a large 
number of multiplications of template values by picture values for each position of the 
template relative to the picture. The process can be speeded up by using Fourier transform 
techniques to compute the cross-correl?tion; however, this approach is more costly in com- 
puter memory requirements, and will not be discussed further here. 

A more quickly computable measure of the mismatch between a template and a picture is 
the integrated (or summed) absolute difference. 
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y>-p(x + u, y + v) 


dxdy 


This measure is zero wh 'n t -ind p are identical (over the area T), Computation of this 
measure involves no multiplications; many fewer arithmetic operations are required than 
in the case of the corre’ation measure (for a quantitative comparison see Bainea and Silver- 
man, 1972). In addition, use of the absolute difference measure makes it possible to take 
advantage of significant shortcuts in the matching process. 


Stoppiiitt Criteria 

A major advantage of the sum of the absolute dii fcrences as a mismatch measure is that 
this measure grows monotonically with the number of points whose absolute differences 
have been added into the sum. This means that if we are looking for a point of best match 
in j given region, we can stop adding up absolute differences in a given position (u, v) as 
soon as the sum has grown larger than the smallest sum sc far obtained in any position 
(u^, ). since we now know that (u, v) ernnot possibly be the position of best match. 

We can do even better if we set up an absolute tltreshold, such that if the mismatch measure 
exceeds & at a given point we will not accept that point as a pouit of match. It is now 
sufficient to add up terms of the sum until 6 is exceeded : when this happens, we can reject 
the given point and move to ^he next point. 

By using both mismatch comparison and absolute thresholding, the process of finding a 
best match (below .he mismatch threshold) can be speeded up considerably. For points 
where the match is poor, the threshold should be exceeded rapidly so that oniy a fraction 
of me terms in the sum need be computed for such points. Even for points where the match 
is good, we need not always compute all the terms of the sum, since we can slop as soon as 
we exceed the smallest previously obtained sum. 


Optimum Matching Sequences 

J( has been pointed out by Nagel and Rosenfeld (1972) that the expected amount of time 
required to exceed a mismatch threshold can be reduced by matching template pixels against 
the corresponding image pixels in a preselected order. To »Hustrate this idea, let us consider 
;i simple example. Suppose that the template and image contain only three gray levels, e.g., 
0, I , and 2, and that, in the given r^ion of tne image, these levels occur with relative fre^ 
quencies of 1 /2, 1 0, and 1 /6, respectively. If • mplate is not at or near i. * correct 
match position, we can assume that the :mage el that coincides with a given point 

of the template is uncorrelated with the template jZky level at that point. Thus, we can 
compute the expected absolute gray level difference between template and image at a point 
as follows: 
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Template Probability of 


Gray Level 


Difrerence 



0 


-1 

0 

1 



1/6 

1/3 

1/3 

0 

0 

1 

0 

1/6 

1/3 

1/2 

0 


0 

0 

1/6 

1/3 

•/2 


Probability of Absolute 
Difference 

Expected Absolute 
Difference 

0 1 



1/2 1/3 

1/6 

2^3 

1/3 2/3 

0 

213 

1/6 1/3 

1/2 

4/3 


In other words, the expected absolute difference between template and image at a template 
point having gray level 0 or 1 is 2/3; but. at a template point of gray level 2, the expected 
absolute difference is 4/3. Thus we can expect, on the average, a faster rate of growth of 
the sum of absolute differences if we compute the difference first for template points that 
have gray level 2. Thus, the mismatch threshold (or the lowest level of mismatch previously 
foun"*) is likely to be exceeded sooner if we compare template points with image points in 
a special order-namely. first using template points of gray level 2 -rather than using the 
template points in an arbitrary order* 

( meraiizing the example just given, it is easily seen that the mismatch threshold (absolute 
or relative) will be exceeded faster on the average if the template points are compared with 
the corresponding image points in decreasing order of expected absolute gray level differ- 
ence. The degree of speed-up achievable in this way depends, of course, on the probability 
distribution of gray levels in ;he given region of the image: if all gray levels are equally 
likely, no speedup is possible by this method. 

Tlie method just described should be even more effective when we are matching midti- 
spectral, rather than gray scale, images. This is because a probability distribution of multi- 
spectral vectc’* values should be even more nonuniform than a distribution of gray levels 
would be. Thus it should be possible to select template points whose expected absolute 
differences from the corresponding image points should be oMite large. Here the dift<^r- 
ences must, of course, be defined in vector icrms, e.g., as sums cf absolute differences of 
the individual vector components. 

The speedup obtainable by this method will not be the same for all landmarks. But the 
best results should be obtained for landmarks in whose vicinity the distribution of gray 
levels, or multispectral values, is higltly nonuniform. Given a set of available landmarks, 
one can determine, for each of them, how much speedup can be expected, and one can 
select preferred landmarks for which the expected speedup is greatest. 


* It may take more time to locate th». Jesired template points than it would if we tested them, perhaps in a simple 
raster sequence; but Nagel and Rosenfeld (1972) have shown that their method is faster than raster matching in 
spite of this. 
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The position of best match, determined by any method, may not be sharply defined, since 
there may be a set of nearby positions for which the match is nearly as good. It lias been 
known for at least 20 years that match sharpness can be improved by differentiating or 
by high-pass Filtering the image and the template before matching them. In other words, 
outlines of regions give sharper matches than do solid regions.* It would be of interest to 
investigate whether the matching speedup process described above is improved by using 
differentiated or outlined images and templates. This will certainly be the case if the 
differentiated image has a less uniform distribution of gray levels (or spectral values) than 
the original image, as we would normally expect. 

ITERATIVE RESAMPLING 

Match positions can be determined to within closer than a pixel by resampling the image 
(at a grid of points that have fractional coordinates, :n terms of the original sampling grid) 
and assigning gray levels to the new pixels by interpolation. The match at such an inter- 
polated position may be better than the match at any of the original (intercoordinate) 
position. The investigation of interpolated match positions is probably best done as a fine 
tuning of a coarse match position found without use of interpolation. The interpolated 
images will tend to have more uniform gray level or spectral distributions than the original 
image, so that the matching speedup process will probably not be as effective during this 
fine tuning stage as at the coarse search stage. 

A refinement allowing more accurate (part pixel) registration, when such accuracy is 
warranted, may be described as follows: 

Let F (x, y) be the cross-correlation function in terms of picture displacements (x, y) from 
the assumed best match point, x = 0, y = 0. Fit a quadratic to F (x, y) of the form 

f y) “ ax^ + 2h + by^ + 2gx + 2fy 

We may determine the coefficients from the values of F at a few points. At x = n, y = 0, 
for example; 


Similarly, 


F(n, 0) = an^ + 2gn 


F (-n, 0) - an^ - 2gn 


hence 


a = { F (n, 0) + F (-n, 0)]/2n^ 
g = ( F (n, 0) - F (-n, 0)] /4n 


* It may be of further advantage to use a thresholded outline image so as to reduce the effects of overall gray level 
differences due to seasonal changes, and so forth. 
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We can also find 


b = [ F (0, n) + F (0. -n)] 2n* 
r = I F(0, n)-F(0,-n)l/4n 
h = [ F(n, n) - F(n. 0) - F (0, n)l/2n^ 

Differentiating. F (x, y) is a maximum when* 

3^0 * hy„ + g = 0 

and 

so that, 

x„= (hf-bgV(ab-h^) 
y^= (gh -af)/(ab -h^) 

Then shift origin to (x^. y^) and repeat with n/2 area size. Since only 5 values of F need to 
be obtained per iteration, the speed of such a procedure may prove to be adequate. 

As an example in one dimension, consider the following template: 

3 19 13 

and the image line gray values: 

6 6 4 8 5 6 6 

Take (F (x)) = ax^ + 2gx. 

(3X 6) + (l X 6)t(9X 4)-t-(l X 8) + (3X 5) 

■v/(3^ + 1 ^ + 9^ + 1 ^ + 3^) (6^ + 6^ + 4^ + 8^ + 5^) 

= 0.6208 = a-2g 

(3 X 4) + (l X 8) + (9X 5) + (l X 6)+(3 X 6) 

F(l) = — - — — . - - 

>/( 3 ^ + 1 ^ + 9 ^ + 1 ^ + 3 ^) ( 4 ^ + 8 ^ + 5 ^ + 6 ^ + 6 ^) 

= 0.6656 = a + 2g 
a = 1/2(F(1) + F(-1)) = 0.6432 
g = 1/4(F(1)- F(-l)) = 0.0112 


After determining coefficients, we require that F be positive definite, i.e., ab > h^. to guarantee a true maximum of F. 
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Differentiating. F (x) is a maximum when ax^ + g = 0, i.e., for x^ = -g/a = -0.01 74 we see; 


Gray Scale 

Old Coordinate 

New Coordinate 

6 

-3 

-2.9826 

6 


-1.9826 

4 

-1 

-0.9826 

8 

0 

0.0174 

5 

1 

1.0174 

6 


2.0174 

6 

3 

3.0174 

In this new coordinate system we will need F(-0.5), hence 
the gray level values at X = -2.5. -1 .5. -0.5, 0.5. 1 .5, 2.5. 

by interpolation we will need 


If linear interpolation is used, then 


(G(Xj)-G(x,)) 

G(x)= G(x )+ ; — (x-x ) 

(Xj-X,) 

thus. 


G(-2.5)= 6.0000 
G(-I.5)= 5.0348 
G(-0.5)= 5.9304 
G ( 0.5) = 6.5522 


and so forth. 

Then a new iteration is undertaken and the process continues until the origin shift indicated 
is satisfactorily small. 

OTHER SPEEDUPS 

A possibility for still further improvement in matching efficiency, not discussed here up to 
now. is to identify locations in the image where the match can be expected to be good. This 
is typically done by making some simple local measurements on the image, and finding 
positions where the values of these measurements are close to their values for the template. 

In the present application, the positional uncertainty of the landmarks is not expected to 
be very large, so that it seems hkely that much will be gained by the use of this approach. 

There is. however, an important possibility for reducing the number of match positions that 
need be tried, once a match has been found f jr the first landmark. If this match is correct, 
the positional certainty of the remaining landmarks has now been reduced, since there 
are now fewer degrees of freedom. This is true even if the match is regarded as only approxi- 
mate. Thus, after the first match is found, the next match can be searched for over a 
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ivlaii\ol\ smaller region. Alter enough matches have been toiind. additional matches liave 
vcr\ little positional uncertainty, so tliat we can search for them over very small regions. 

Of course, if the subsequent matches are not found in the expected places, they must be 
sought for over wider regions, and tlie original matches must be reevaluated. Tints a feed- 
back process can be used to zero in on a best combination of match positions. 

TECHNICAL DISCUSSION 

The specific application to which tnis rcscarcli is directed is the navigation of the Synchronous 
Meteorological Satellite (SMS I. In umc ii is lioped that the processing of recognized land- 
marks for navigation (both for location error correction and relocation) can take place in 
quasi-real-time, which in this context is approximately 5 minutes. Before this can come 
about, there is a need to determine landmark processing algorithms that are rapid, accurate, 
reliable, and economical in storage. 

It is recognized that several operational questions influence tlie utility of the algorithm 
development results. In particular, eventual implementation may be effected on a mini- 
computer with a hardware dot product. The availability of such a feature would have 
some significance to the efficiency of algorithms in the multispectral domain. 

Also, error budgeting will play a central role. It may not be useful, for example, to process 
the image data with extreme accuracy since the errors in the landmark surveys may well 
predominate (current landmark accuracies are believed to be about 1/iOOth of a degree in 
latitude and longitude, and may have a systematic bias). Hence it may be useful to con- 
centrate on using the same landmarks repetitively for relative image -to-i mage geometric 
transformations rather than absolute geodetic location determination. 

Thirdly, cloudcover distribution probability functions will affect the ability to use the land- 
marks. The geometric relationships require a well distributed set of control points if orien- 
tation and orbit parameters are to be determined accurately. Inadequate or ill-distributed 
control point data give rise to ill-conditioned matrices and the calculations become unreliable. 
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CLOSED FORM SATELUTE THEORY WITH 
EQUINOCTIAL ELEMENTS 


R. Broucke 

University of Texas 
Austin. Texas 


A computer software system has been developed to generate the closed-form literal first- 
order perturbations due to any harmonic in the potential. In a first approach, classical 
elements and the true anomaly are used. In another approach, equinoctial elements and the 
true longitude are the basic variables. The solution in equinoctial elements does not have 
any zero eccentricity or zero ir 'ination singularities. 

In the case of tesseral and zonal harmonics, the rotation of the central body is neglected. 
The expansion of the potential is done with simple recurrence formulas. 
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FORMULATION OF AN ARBITRARY GEOPOTENTIAL TERM 
(TESSERAL) IN EQUINOCTIAL VARIABLES 


P, Cefola 

Computer Sciences Corporation 
Silver Spring, Maryland 


Previously, general formulas for the averaged disturbing potential were obtained ( AlAA pre- 
print 75-9) in equinoctial elements for the zonal and third-body harmonics, in addition, 
methods were given for recursive computation of these potentials. The current paper extends 
the applicability of this model by giving an explicit expression in equinoctial elements for 
the arbitrary geopotential term: 



P (sin 6 ) (C cos m X + S _sin m X] 

nm ‘ nm nm * 


Expressions for the spherical harmonics P^ ^ (sin 0) cos mX and P^ ^ (sin <t>) sin mX are 
obtained in terms of Jacobi polynomials with the argument (1 - p^ -q^ )/(l + p^ + q^ ) 

(p and q are equinoctial elements), polynomial functions which are straightforward gen- 
eralizations of the and polynomials appearing in the zonal potential, the true longi- 
tude, and the Greenwich sidereal time, 6 is fixed during the averaging period. However, 
an expansion of the true longitude in terms of the mean longitude and the equinoctial 
elements k and h would allow the consideration of resonant cases. Finally, consideration 
is given to recursive computation of the averaged potential for tesserals. 
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AN ANALYTIC METHOD TO ACCOUNT FOR ATMOSPHERIC 
DRAG WITH AN ANALYTIC SATELLITE THEORY 


N. Bonavito and R. Gordon 

Goddard Space Flight Center 
Greenbelt, Maryland 


The motion of an artificial earth satellite in the presence of air drag and the eartlTs gravita- 
tional potential is considered. In contrast to the classical methods of numerical integration, 
this approach presents a quadrature algorithm employing analytical expressions for the 
variation of orbital elements produced by the air drag. These expressions are w'ell-defined 
over expanded subintervals of the solution, and produce accurate agreement with profiles 
of tabular density. This procedure then allows a flexibiliiy in the selection of end points 
of the subintervals, which in turn ensures a minimum error bound on the required analytical 
function. 

In this method the effect of oblateness is accounted for by either the Vinti spheroidal theory, 
the Brouwer orbit theory, or the Brouwer-Lydanne theory. The changes due to atmospheric 
resistance for a nonrotating sphere are accounted for by the solutions of the variational 
equations, evaluated with the appropriate theory. 
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NUMERICAL IMPLICATIONS OF STABILIZATION 
BY THE USE OF INTEGRALS 


P^ul R. Beaudef 

Computer Sciences Corporation 
Silver Spring, Maryland 


The subject matter of this paper is some numerical experiences that we iiave had inx o!' * 
some of the celebrated notions of dynamic stabilization. Figure 1 demonstrates Ljapunov 
stabilization, which is an analytic notion, and die example given is one of circular motion 
under an arbitrary attractive central force field. There is some particle moving around in 
a circle at some distance from an attractive center. The attractive force is designated by T(r). 
the velocity of the spacecraft is v, and the centripetal acceleration is v^ /r. By equating the 
centripetal acceleration to the attractive force, we can get an expression for the angular 
frequency that is the square root of the attractive force over r(y^f/r). 

Notions of stabilization involve questions associated with what happens to this motion under 
slight perturbation of initial conditions. Let us start out with some satellite or particle 
moving in a circular orbit at a radius of a; there will be a certain angular rate associated with 
that motion. If we cause a slight change to occur in the initial conditions so that the radius 
of the orbit is no longer a but is a + Aa, then it may happen that the angular rate may differ 
from that of the original orbit. These initial conditions should be selected in such a way that 
the motion will still be circular in this example. 

As a result of the possible different rates in the angular frequency, the mean anomalies (0\) 
between both satellites will be different and their difference will increase linearly in time. 
Such motion is dynamically unstable because the motion of the initial problem and that of 
the perturbed problem (that is, the problem with slightly perturbed initial conditions) will 
deviate arbitrarily; the deviation in the mean anomaly will be as great as desired if a 
sufficient length of time is allowed to elapse: 


for any 



(0-^2 (t)>6 


for sufficiently large t. 
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V = VELOriTY 


2 



v^/r = 

centr:petal 

ACCELERATION 


-f(r) = 
ATTRACTIVE 
ACCELERATION 



Figure 1. Ljapunov stability, circular motion under arbitrary attractive central force. 

The single exception to this rule is the case where the angular frequencies are in fact the 
same for all possible radii. This implies that the force taw should be given by -w^r, which 
is the simple harmonic oscillator; 

cj| = cOj = CO -*■ f (r) = - r. 

Since the orbit (Keplerian) problem involves the force law, which is inversely proportional 
to the square of the field, this i..iplies that the Keplerian problem is Ljapunov unstable: 

f(r)= 

r^ 

A dynamic problem is going to be either Ljapunov stable or unstable This is a physical 
concept and depends on the particular problem. The notion of stabilizing an unstable prob- 
lem must then involve a change in the problem itself. We must look for another problem 
that happens to have the same solution as the original problem, and how this can be done 
will be explained later in the presentation. 

The reason for looking at analytic concepts of stability is to try to improve the accuracies 
associated with numerical integration. Two basic approaches have been examined: With 
the two-body motion, one possibility (A) is to find a different formulation of the equations 
of motion, which are dynamically stable (that is, Ljapunov stable). Such examples are 
given by Baumgarte or the simple harmonic oscillator Kustaanheimo Stiefel (KS) theory. 
For such formulations, the frequency of the motion— the energy— is an a priori constant. A 
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second procedure (B) is to add a constraint (for example, an energy constraint) in addition 
to the equations of motion. 

H (X,X) = - Pq = constant 

If there were such a constraint, then, when initial conditions were perturbed, those initial 
conditions would have to satisfy the constraints. Variations in the initial conditions that 
change the energy (in the case of frequency constraint) would not be permissible alterations 
in the initial conditions. 

So it happens that, for the tv/o-body problem, concepts A and B are analytically, but not 
numerically, equivalent, as will be shown later. 

With perturbed motion, the Baumgartc and KS theories are not Ljapunov stable. In fact, 
they are unstable, but the degree to which this instability occurs is of the order of the 
perturbations and numerically this does not present any difficulty (concept A). On the 
other hand, employing an energy constraint (B), the negadve energy may or may not be 
a constant. If we had a problem, the Keplerian energy would be a constant, yet the 
Baumgarte or the KS theories would still be unstable. So, for perturbed motions, concepts 
A and B are not equivalent, and the question arises as to whether it is better for numerical 
integration to look for Ljapunov stable or nearly Ljapunov stable equations of motion or to 
apply some kind of energy constraint. All of the evidence indicates that slightly better 
answers are obtained, at least for nearly circul r motions, by using energy constraints. 

We now take a look at some of the methods of dynamic stabilization, by which is meant 
either formulations that are Ljapunov stable or nearly Ljapunov stable and/or applying an 
energy constraint. The crux of all the stabilization procedures that exist are presently being 
investigated: The Baumgarte approach; a recent approach that was proposed by Stiefel, 
with which we have not yet had any numerical experiences; the Baumgarte-Stiefel stabiliza- 
tion procedure, which is an energy-constraint-type formulation; nd a formulation of 
applying an energy constraint directly, which is due to Nacozy and will be discussed later. 

The Baumgarte procedure involves a Ljapunov-stable system of equations for the two-body 
problem and is based on a method of Poincare. If we start with a Hamiltonian, which is a 
function of the position and conjugate momentum (velocities), then the Keplena.n equations 
of motion are given by the following canonical equations: 

8H dH 

H(X,P)^X= and P = - 

dP ax 

In numerical integration, it is very often desirable, especially for eccentric motion, to make 
a transformation from time as an independent variable to some other independent variable, 
an S variable. This S variable could be the eccentric anomaly for certain eases or a mean 
anomaly in another case, but it does not matter. The type of equation which will relate time 
with the new independent variable is given by dt/ds as some function and, in general, that 
function may be a function of both position and momentum: 
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dt 

lls 


Poincare defined a new Hamiltonian, which is related the old Hamiltonian by the func- 
tion (dt/ds) times the old Hamiltonian plus a constant, and it turns out that, wiih this new 
Hamiltonian, the equations of motion for both the space variables and the ti.. e variable 
are canonical, as shown below: 


, dH m o/i 

X' - =fi f — (H^PJ 

dP dP dP ^ 


dH dH dp 


dX 


dX dX 


dH dp 

ir 




dt 


= 0 


= constant sudi that (H + P^^J = 0 


Tlieie are three equations here that are the correct equations for dynamics, except for cer- 
tain terms w^ich are called control terms. The equation for P^, which is the momentum 
conjugate to the time, defines The S-derivative of P^ should vanish, and this implies that 
P^j is a constant. If that constant is chosen to be equal to the negative of the energy or the 
value of the Hamiltonian, then the control terms would be numerically equal to zero, but 
dynamically they are not. Dynamically, they are some function of r and the constant 
energy. So this ,new canonical system of equations is really an entirely diffeic:)t problem 
that happens to have the same solution as the old problem. 

Baumgarle noticed that for the function p = r the equations of motion with respect to the 
new, indepeiident variable S (that is, the spatial equation of motion w th respect to the new, 
independent variable S) are dynamically stable in the Liapunov sense and, in fact, can be 
transformed to the KS simple harmonic oscillator equations. Tlie time equation. T' = r, 
is not dyiiamically stable, but tfierc are ways of getting around that problem as will be 
shown later in the p iper. 

:n though tiK control term, H + , analytically equal to zero from the analytic solution, 

when these cqiMtions of motion uit numerically integrated, the control term may develop 
some error. We will then want io knovV how this error is going to beb; ve dynamically. 
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To examine this formulation in the time domain, we use the following equation: 


- ^ X 

X=-— X+— [HrPJ 


CONTROL TERM 


We can derive an equation of motion for the control term by multiplying (dot product) 
this equation by the velocity and doing some manipulation. The general solution of this 
equation has the control term as a constant multiplied by r. For circular motions this 
wovld basically be a constant. This solution tells us that if we develop some numerical 
(nonzero) errors in the control terms, then those errors are going to persist. Hence, from 
a numerical point of view, the orbital state will not be on the correct energy surface. 

This is different from our analytical approach, in which we saw that, for the two-body 
problem, Baumgarte stabilization (or Ljapunov stabilizatic*') and applying an energy con- 
straint are, in a way, equivalent. From the numerical point of view we see that this is not 
true, because in numerical integration we create errors ^nd those errors persist in the 
BaumgLite approach, so the state is not forced back onto the correct energy surface. Thus, 
applying the Ljapunov stable system of equations does not ensure that the state is going io 
be on the correct energy surface, 

i here are otiier methods of stabilization. The Stiefel approach involves multiplying the 
forcing terms of the differential equations by the a priori energy constant divided by the 
Hamiltonian: 


Liiven 

X = f(Zo 


use 



Now this parentheses is ideally equal to 1 for the analytic solution. If we raise the exponent 
to the power, numerically it does nothing, but from a stability point of view, it does. 
This Stiefel system of equations is dynamically stable in the Ljapunov sense. 

The baumgarte-Stiefel control term is similar to the Baumgarte control term, except that 
the muUipiie is proportional to the velocity vector over the square of the velocity The 
multiplier is sometimes called a dissipative term, because,, by looking at the eauation of 
motion that governs the control term, it ca *’ seen that it has a solution which decays in 

time: 
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leads to 


K^X 
X = 


— IH + Pol 


_d 

dt 


'H + Pol 


-^[H + Pol 


[H + P^jl = constant 


Hence, the Baumgarte-Stiefel control term forces the state back onto the correct energy 
surface, because the control term is forced to decay to zero. 

In the Nacozy approach, at each step of the integration, the control term is forced to zero 
by a virtual displacement made in the state, such that the new energy after such a displace- 
ment should be the correct energy. 

6X = - F = acceleration. 


It happens that the Nacozy procedure is not compatible with multistep integration proresjts 
because of the di*^continuity that occurs in the state variables. The numerical results 
single-step numerical integration methods apparently work; with the multistep metho/i^. 
there are problems, as will be shown later. 


We now turn to tht p*»rturbed problem where we have nonconservative perturbation, such 
as drag or solar radiation pressures, conservative forces, which are derived from the gracing 
of some potential, and control terms: 


X = 


^ ^ - 

X + P “ W + [control term) 


The^e are two types of integral onstraints to discuss: The Keplerian eneigy is ? near integral 
of ‘he motion and has the equation: 


dt 




• P 


There is also a term not ino . ■ here, which is the partial of the potential with respect to 
time. the absence of perturbation, this 's zero, so it is a near constant of the motion. 
Instead of having a that is e const.^nt, *ve have to numerically integrate it. The following 
is a comparison of that integration of with the computation of the Hamiltonian in a 
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conservative system. The sum of H aru: is ideally zero and will be used as a control term: 



Tlie other type of integral constraint is a Jacobi integral, which is what the energy looks 
like in the rotated coordinate system of the earth. It is an exact constant of the motion 
tor the TuH potential including tessera! harmonics: it is given by the Hamiltonian for the 
inertial energy plus the angular velocity of the earth dotted into the angular momentum. 
When V is time-dependent due to the earth’s rotation (tesseral harmonics): 

— [H + w -(RX R)J = + (X-Z5 X R)- P =0IFP=0. 
dt 

A control term based on this integral of motion can be applied, designated in later equations 
by R, indicating the energy in the rotating coordinate system: 

[H + (o . (R X R) + P^] 


where 


dt 


-(X-ojX R)- P 


None of these equations are integrated in the time domain. They are integrated usi g the 
eccentric anomaly as an independent variable and, in addition, the time equation must be 
integrated. 

Baumgarte dynamic stabilization also requires a time equation that is dynamically stable. 
There are five options for the time equation. The normal one would be t' = r, the defining 
relationship between t and s. 

If we differentiate this equation with respect to s, we get r" - r'/K. The reason we consider 
this is that second-order systems of equations are often easiei to integrate using class 2 
methods. 

The next option is the ordinary time element option: 


t = T- 


Kr' 


2Po 
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The equations of motion are essentially a constant term plus perturbation: 


K 



IR • (P - VV) - 2V| 



This is a constant when the energy is a constant. If the energy has to be integrated, 
then some errors in P^ might be expected to crop up: if those errors are significant, then 
there can be secular error type terms that will grow. There will then be timing errors 
associated with the integration of 

This analysis leads to a third system of equations, in which a new eccentric anomaly is 
related to the old eccentric anomaly by the following combination of variables: 

dn K d K d 

ds 2Pjj ds 2P^ dT 7 

Here the energy P^ is in the denominator. We can traiisform all of the equations of motion 
by using this relationship. 

The time is then related to a new time element, t*, by adding to it a modified eccentric 
anomaly; as the equations of motion for t* involve only perturbations, in the absence of 
such perturbations, the anomaly is zero: 

Kr' 


dr* r 1 ^ 

— = — [R*(P -W)-2V) - 

df? K ^ Po 

The fifth possibility for deriving the time as a function of s is to use a third-order differential 
equation as shown here: 

t"’ + 2P,,t'=K 

We have not had any numerical experiences doing this, but this equation is dynamically 
stable, as are the time element equations. 

Tables I through 7 give the numerical errors of the GEOS-B test orbit from computer sim- 
ulations. The errors ,^re given in units of 10 * kilometers (millimeters) after integrating the 
GEOS-B orbit over 50 revolutions. The perturbations are a 1 5-by-l 5 geopotential field, 
and the results are given over a range of 60 to 100 integration steps per orbit. In the first 
column, an 1 indicates the integration of the inertial energy, and an R indicates the integra- 
tion of the rotational or body-fixed energy. The second column provides the particular 
time equation used to get the time from the eccentric anomaly. A zero indicates the t" 
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time equation, 1 corresponds to the t' equation, 2 is the ordinary time element equation, 
and 3 is the modified time element equation. The last four columns give the errors for the 
Baumgarte dynamic stabilization procedure (B); the Baumgarte-Stiefel control term, which 
forces the state cack onto the correct energy surface (BS); the experiences using the 
Nacozy approach (N); and the ordinary time-regularized formulation, which does not have 
any control term and is dynamically unstable (TR), respectively. 

Related to these tables are the experiences associated with integrating Cowell’s equations 
of motion using time as an independent variable. Since the motion is nearly circular, 
equal step sizes in time will correspond very closely to equal step sizes in eccentric anomaly. 

1 the last four columns of the tables, the three numbers in each entry correspond to intrack 
atial errors, crosstrack spatial errors, and a timing error (which is also an intrack error), 
respectively. There are three errors here because we are integrating spatial equations and a 
time equation with respect to eccentric anomalies. The final time at the end of the run is 
one in which the eccentric anomaly has reached a final fixed value at the end of 50 revolu- 
tions. The results of these runs are as follows: The time element with the Baumgarte con- 
trol terms works with the timing error and the intrack error is still somewhat large compared 
to the crosstrack error. This result might be expected, because the state is no longer on the 
exact energy surface but persists off the energy surface. When a modified time element is 
applied, the timing error almost completely disappears, but only at the expense of intrack 
errors in the spatial equation. So all the modified time element has done for us is to shift 
the timing errors to spatial errors. The Nacozy process has really done nothing for us 
either, and the reason for that is the discontinuities in the state. Taking finite differences 
just adds errors in a multistep integration process. If we had a single-step integration process, 
these problems would disappear, as other experiences have indicated. The time-regularized 
processes are unstable and about two-thirds of the intrack error is in the timing and about 
one-third is in ihe spatial equations. This is in agreement with analytic work that has been 
done by Baumgarte. 

It is the Baumgarte-Stiefel control term that gives the best answers. If, for example, the 
Jacobi integral is used, we get very little timing error. The spatial errors are still a little bit 
larger, but this control term gives us the best overall results. An important note is that the 
Baumgarte-Stiefel control term does not require a time element. In fact, the best results 
were achieved when the time element was not used and when the T" equations or the T' 
equation were integrated. This is very important. Basically, the timing error results from 
consistent errors in r; T' = r. If there are consistent errors in r, timing errors will develop 
in integrating that system, the double integral problem. 

By applying the energy constraint, the radius r is adjusted so that consistent errors in the 
time do not occur, and this has led to the results presented here. It is recommended that 
an energy constraint can be applied without worrying about the timing problem at all. 
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Table I 

Numerical Errors of the GEOS-B Test Orbit 
for the BDSP- 60 Steps Per Orbit 


ENERGY 

ITELEM 

B 

BS 

N 

TR 

1 

0 

UNSTABLE 

-127.07899 

28.97938 

561.96729 

-5802.32261 

84.70741 

12330.44769 

-5326.43758 
13-81227 
12123 72627 

1 

1 

UNSTABLE 

-127.07197 

28.97945 

282.55283 

-5802.18569 

84.70673 

6753.60074 

-5326,43611 

13.81335 

11831.B5920 

' 

2 

UNSTABLE 

-127.11421 

28.96996 

244.29093 

-5803.50299 

84.74010 

1405.50968 

-5327.36544 

13.81054 

665.25064 

1 

3 

UNSTABLE 

-341.064 

28.773919 

42.056 

5693.628 

87.159853 

1136.962 

-5244.249 

6.950228 

438.872 

R 

0 

UNSTABLE 

- 22.85710 
29.93991 
331.36557 

UNSTABLE 

SAME AS 10 

f' 

1 

UNSTABLE 

- 22.84957 
29.93998 
50.23874 

UNSTABLE 

SAME ASH 

R 

2 

UNSTABLE 

- 22.84268 
29.95027 
12.51864 

UNSTABLE 

-5327.29464 

13.81072 

1799.09322 

R 

3 

UNSTABLE 

- 22.84268 
29.95027 
12.41565 

UNSTABLE 

-5327.29408 

13.81074 

1798.98859 


NOTES: 

In the last four columns, the three r imbers in each entry (given in units of 10*^ km) 
correspond to intrack spatial errors, crosstrack spatial errors, and a timing error (which is 
also an intrack error), respectively. 

By integrating Cowell’s equations of motion with time as an ind mdent variable, the intrack 
spatial error was 3401 .4027, and the crosstrack spatial error was 26.77772. 
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Table 2 

Numerical Errors ol the GEOS-B Test Orbit 
for the BOSP. 65 Steps Per Orbit 


ENERGY 

ITELEM 1 B 

BS 

N 

TR 

1 

0 

-114.19500 

2.53846 

1828.84252 

-115.66616 
14 78646 
33948653 

-2006.31579 

4574424 

4266.29205 

. 

-1814 58633 
1.25166 
4120 60271 

1 

1 

-114.19606 

2.53844 

1867.57285 

-115.66429 

14.78647 

261.06752 

-2006.27240 

45.74470 

2416.86411 

-1814.58677 
1 .25165 

4031.98843 

- -- 

1 

2 

-114.18013 

2.53900 

204.35943 

-1 15.69274 
14.79099 
259 28441 

-2007.21942 

45.76519 

648.56484 

-1815.31460 
1.25176 
399 31962 

1 

3 

-337.19 

3.09228 

11.296 

-338 474 
14 558205 
43.667 

2117.405 
46 343831 
409.824 

1936.262 

3.169659 

175.349 

R 

0 

-101.88979 

2.52544 

1825.45528 

- 29 81388 
15.44724 
148.96701 

UNSTABLE 

SAME AS 10 

R 

1 

-101 89076 
2.52543 
185842729 

- 29.81174 
15.44727 
7.00168 

. j 

UNSTABLE 

SAME AS 11 

R 

2 

-101.87554 

2.52593 

307.78874 

1 

- 29.80779 j 
15.45175 
65 45256 | 

UNSTABLE 

i 

-1815.25454 

1.25175 

667.86354 

R 

] 

3 

-101.87549 

2.52593 

307.SS044 

29.80778 ' 

15.45176 I 
65.31316 

i 

UNSTABLE 

-1815.25555 

1.25173 

667.72535 


NOTES: 

In the last four columns, the three numh«*N in each entry (given in units of 1 0'^ km) 
correspond to intrack spatial errors, c track spatial errors, and a timing error (which is 
also an intrack error), respectively. 

By integrating Cowell’s equations of motion with time as an independent variable, the intrack 
spatial error was 3401.4027, and the crosstrack spatial error was 26.77772. 
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Table 3 

Numerical Errors of the GEOS-B Test Orbit 
for the BDSP, 70 Steps Per Orbit 


ENERGY 

ITELEM 

B 

BS 

N 

TR 

1 

0 

- 56.39602 

3.11021 

- 83.00906 

- 55.74016 
3.72737 
110.25181 

89.17622 
9.39301 
- 190.00997 

59.21339 
3.54422 
- 144.20532 

1 


- 56.39740 

3.11020 

- 32.41277 

- 55.74068 
3-72737 
126.90710 

89.17136 
9.39302 
- 39.18658 

59.21225 
3.54421 
- 130.56524 

1 

2 

- 56.41604 
3.10982 
106.19045 

- 55.75580 
3.72860 
136.67825 

88.64037 

9.40287 

109.40972 

58.77161 

3.54320 

129.18683 

1 

3 

-172.115 

2.960583 

2.72 

-171.421 

3.569666 

25.41 

1.232 

8.767146 

7.63 

31 .457 
3.383486 
15.855 

R 

0 

- 18.50051 
3-14627 
-127.93949 

- 22.49304 
3.90129 
36.38436 

UNSTABLE 

SAME AS 10 

R 

1 

- 18.50192 

3.14625 

- 77.63998 

- 22-49356 
3.90129 
52.98888 

UNSTABLE 

SAME AS t1 

R 

! 

2 

- 18.518117 
3-14591 
15.28445 

- 22.49164 
3.90259 
60.43499 

UNSTABLE 

58.81103 

3.54332 

35.63308 

R ! 

3 

- 18.51815 1 

3.14591 
15.35058 

- 22.49178 
3.90258 
60.50244 

UNSTABLE 

58 80972 
3.54331 
35.70119 


NOTES: 

In the last four columns, the three numbers in each entry (given in units of 10 ^ km) 
correspond to intrack spatial errors, crosstrack spatial errors, and a timing error (wliich is 
also an intrack error), respectively. 

By integrating CowelPs equations of motion with time as an independent variable, the intrack 
spatial error was 3401.4027, and the crosstrack spatial error was 26.77772. 
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Table 4 

Numerical Errors of the GEOS-B Test Orbit 
for the BDSP, 75 Steps Per Orbit 


ENERGY 

ITELEM 

B 

8S 

N 

TR 

1 

0 

- 19.$il49 
Y. 84296 
-545.21889 

- 19-12725 
.99247 
10.04084 

631.43636 

5.75062 

-1342.56385 

534.84907 

2.43962 

-1222.06456 

1 

1 

- 19.62227 
1.84295 
-515.99802 

- 19.12819 
99247 
44.14811 

631.41960 
5-75071 
- 698.77660 

534.84874 

2.43962 

-1187.46233 

1 

2 

- 19.64372 
1.84251 
41.23325 

- 19.13503 
.99227 
53.5391 1 

631.16381 
5.74633 
- 77.89860 

534.62285 

2.43881 

9.64151 

1 

3 

64.153 

1.592328 

.078 

1 

63.651 

5224b 

12.097 

588.847 
5.790954 
- 119.748 

492.282 

1.808203 

31.944 

R 

0 

9.20689 i 
1.8C947 
-581.840 

- 13.39047 

1.023638 

- 2.75868 

UNSTABLE 

SAME AS 10 

R 

1 

9.20615 1 

1.86948 1 

-551.67891 ! 

- 13.39146 
1.02364 
31.42207 1 

. 1 

UNSTABLE 

SAME AS 11 

R 

2 

9.18710 
1.86908 
- 64.66871 

1 

- 13.39082 1 

1 .02341 < 

39.46957 

UNSTABLE 

i 

534.64250 
2.43889 
- 142.20691 

rt 

3 

9.18710 
1.86908 
- 64.61750 

- 13.39086 
1.02342 
39.51002 

UNSTABLE 

i 

1 

534.64266 
2.43889 
- 142.15558 


NOTES: 

In the last four columns, the three numbers in each entry (given in units of 10 * km) 
correspond to intrack spatial errors, crosstrack spatial errors, and a timing error (which is 
also an intrack error), respectively. 

integrating Cowell’s equations of motion with time as an independent variable, the intrack 
spatial error was 3401 .4027, and the crosstrack spatial error was 26.77772, 
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Table 5 

Numerical Errors of the GEOS-B Test Orbit 
for the BDSP, 80 Steps Per Orbit 


EN6RGY 

ITELEM 

CD 

— 

6S 

— 

N 

TR 

1 

0 

- 3.88340 

.85629 
-459.15305 

- 3.74365 
1.48006 

- 15.88953 

547.01168 

7.57810 

-1161.99962 

462.17194 

1.23206 

-1052.62631 

1 

1 

- 3.88368 
.85629 
-448.19016 

- 3.74432 
1.48006 
9.26784 

546.99724 
7.57919 
- 619 64184 

462.17242 
1 .23206 
-1026.10445 

1 

2 

- 3.89834 

.85601 
10.98153 

- 3.74692 

1.48023 
15.50786 

546.89183 
7.57760 
- 96.19190 

462.07088 
1.23162 
- 21 .79383 

1 

3 

- 16.867 

.699993 

.639 

- 16.724 
1.182493 
4.999 

530.505 
7.653788 
- 105.966 

446.289 

.752743 

32.042 

R 

0 

12.32661 

.86888 

-480.49666 

- 6.75389 
1.54807 

- 9.24393 

UNSTABLE 

SAME AS 10 

R 

1 

12.32631 

.86888 

-468.66261 

- 6.75457 

1.54807 
15.97804 

UNSTABLE 

SAME AS 11 

R 1 

2 

12.31317 
.86863 
- 59.98543 

- 6.75434 

1.54821 
21.53818 

UNSTABLE 

462.07985 
1.23164 
- 134.05535 

R 

3 

12.31314 
.86863 
- 60 08864 

- 6.75426 
1.54821 
21.43552 

UNSTABLE 

462.079696 
1.23164 
- 134.15899 


NOTES: 

In the last four columns, the three numbers in each entry (given in units of 10 * km) 
correspond to intrack spatial errors, crosstrack spatial errors, and a timing error (which is 
also an intrack error), respectively. 

By integrating Cowell’s equations of motion with time as an independent variable, the intrack 
spatial error was 3401.4027, and the crusstrack spatial error was 26.77772. 
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I cihlc (> 

Numerical I rrors ot the (]| OS-B Test Orhit 
tor the BDSP> ^>0 Steps Per Orhit 


ENERGY 

ITELEM 

B 

— 

BS 

— 

N 

TR 


0 

1.94060 

.17462 

-135.07206 

1.88804 
.65661 
- 9.08666 

172.74506 
3.25877 
- 365.23276 

144.39606 
.21146 
- 326 48960 

1 

1 

1.94069 

.17463 

-137.79058 

1.88791 
.65661 
- 3.40434 

172.74175 
3.25882 
- 198.53333 

144.39652 
.21147 
- 319.98778 

1 

2 

1.93719 
.77468 
- 1 .66455 

1.88796 
.65678 
- 1.60664 

172.73643 
3.25888 
- 36.33632 

144 38806 
.21146 
~ 1310179 

1 

3 

1.632 

.212515 

.340 

1.581 

342734 

.362 

170.214 
3.062124 
- 33.904 

142.24 
-192607 
- 10976 

R 

1 

0 

4.84812 

.17485 

-139.27516 

.91961 
.68519 
- 2.84861 

UNSTABLE 

SAME AS to 

1 

R 

1 

4.8481 

.17485 

-141.66259 

- .91983 

.68519 
2.84094 

UNSTABLE 

SAME AS <1 

R 

2 

4.84502 
.17489 
- 18.61224 

91976 

.68634 

4.51218 

UNSTABLE 

144.38850 
.21146 
- 43.66740 

R 

3 

4.845096 
.17489 
- 18.82819 

.91971 
.68534 
- 4.29666 

UNSTABLE 

144.38868 
.21146 
- 4388346 


NOTPS: 

In the last four columns, the three numbers in each entry (given in units of 10'^ km) 
correspond to intrack spatial errors, crosstrack spatial errors, and a timing error (which is 
also an intrack errorf respectively. 

By integrating Cowelfs equations of motion with time as an independent variable, the intrack 
spatial error was 3401 .4027, and the crosstrack spatial error was 26.77772. 
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l ablo 7 

Numoricul I rrors of the (W OS*B l est Orbit 
for the BDSP, 1 00 Steps Per Orbit 


ENERGY 

ITELEM 

B 

BS 

N 

TR 

, 

0 

1.00589 

.99013 

19.86950 

16.67995 



.078112 

.15110 

.65629 

072534 



- 10.73311 

.96997 

- 40.52781 

~ 3609153 

1 

1 

1.00598 

.99015 

19 86833 

16 68023 



.078109 

.15109 

.65628 

.07253 



- 13.79768 

- 1.27553 

- 22.83707 

~ 36 10359 

1 

2 

1.0 .97 

.99046 

- 19.87749 

16.68774 



.078140 

.15116 

.65649 

072546 



.59675 

.92698 

~ 478582 

2 22072 

1 

3 

.775 

.790 

17.38 

14.315 



.241308 

.171367 

.348605 

.236638 



.155 

.180 

3 891 

1.423 

R 

0 

.83340 1 

.24629 

UNSTABLE 

SAME AS to 



.078159 

.15553 





- 10.70273 

68240 



R 

1 

8C356 

.24624 i 

UNSTABLE 

SAME AS i1 



.078162 

.15553 





~ 13.68274 

.38025 



R 

2 

.83348 

.24621 

UNSTABLE 

16 68663 



.07819 

.15559 


.072566 



^ 1.26500 

.71750 


4.77124 

R 

3 

.83348 

.24632 

UNS' ABLE 

16 68570 



.07819 

.15559 


072560 



- 1.31470 

.66760 


4.82004 


NOTl :S: 

In tlie last four columns, tlie tliree numbers in each entry (‘j^iven in units of 10 ki]i) 
correspoiui to intrack spatial errors, crosstrack spatial errors, and a timing error (which is 
also an inlrack error), respectively. 

By integrating (owelTs cHiuations of motion with time as an imlependent variable, the inlrack 
spatial error was 3401 .4027, and the crosstrack spatial error was lh.1'^111. 
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DISCUSSION 

VOICE: In th<. ..lulation of the Bauingartc-Stiefel stabilization procedure, corrections are 
applied in only one direction of a six-dimensional manifold. Wouldn't you expect more than 
errors in time? Is all the instability controlled by just the energy constraints? 

BEAUDiE: Yes. 

VOICE: It seems that you ought to do more than just force your solution back into the 
correct energy surface. 

BE A UDET: That is in fact correct, but when you look at it, you have to ask yourself the 
question: “What direction in this six-dimensional slate manifold corresponds to the unstable 
direction, t!ie direction in which errors are going to consistently grow?" 

It's like that first diagram I gave you (figure 1 ). It is the fact that the freciuency is in error 
that gave rise to an ever-increasing mean anomaly. It is an intrack error. In other words, 
when we integiate, we don't worry at all about crosstrack radial errors. It's intrack errors 
that are going to grow in time. 

It turns out that the uncertainties in the frequency of the motion cause these errors and, 
since the frequency is somewhat related to the energy, at least in the two-body problem, 
applying an energy constraint solves the stability problem. 

VOICE: But there is no constraint that will give a stable manifold, except for very special 
problems, which indicates that most problems should be »:nstdble. 

BEAUDET: In die perturbed problem that is the case and, if we apply perlurbati(m. 
we're going to have an instability associated with the rotating plane ol the orbit. In tfiat 
direction, the manifold space is going to be unstable, even though we might apply a: .nerg\ 
constrai We’re simply fixing up the worst part of the tilings associated with two-body-lype 
motion. 
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IV/C£: What have yon done to take care of the high drag case? 

BEAUDET: I haven’t added anything into the energy to take care of drag. We have integrated 
orbits where drag was a perturbation and in which I knew that drag was not too severe a 
perturbation. We still get stabilization as a result. Of course, we have to integrate the energy 
equation, and drar a ppears on the right-hand side. The question is always associated with how 
accur lely we can ^grate this energy equation. 

VOICE You mean you can stabilize the problem in the case of drag? 

BEA UDET: If the drag is not too severe. If the drag becomes very severe, we’ve come 
across another problem, that the t centric anomaly is not the right independent variable to 
u«e. If we start hitting a big wall associated with drag, we would like to use a different 
independent variable than eccentric anomaly while integrating through that region. It is 
the difficulty associated with finding such an independent variable or using it that has given 
rise to our inconclusive results in high drag cases. 


c:. % 
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SPECIAL PERTURBATIONS USING BACK-CORRECTION 
METHODS OF NUMERICAL INTEGRATION 


Teny Feagin 

University of Tennessee Space Instate 
TuUahoma, Tennessee 


Previous speakers have discussed how to change the analytical framulation of various prob- 
lems in order to introduce stabilizing effects or to imivove the accuracy and thereby in- 
crease the efficiency of numerical integration. The present topic is concerned with how to 
improve the process of numerical int^ration itself, in m'der to increase efficioicy in the 
development of accurate ephemerides for earth satellites. In particular, the subject of this 
discussion is a new class of linear multistep methods for the rumerical integration of ordin- 
ary differential equations. These methods are distinguished from the classical methods in 
that they permit the solution to be corrected at cotain “back” points. That is, in the case 
of satellite computations, the solution is corrected at certain points in the past as the 
int^ration advances in time. Algmithms have beoi developed for the solution of both 
first- and second-order differential equations, althou^ only the second-order case is con- 
sidered here. 

There are two reasons for correcting the solution at back points: First, when a polynomial 
obtained from interpolating evenly spaced data is used for approximating a function at a 
point, the coefficient of the errm- term is smaller when the point is nearer the middle of the 
range of data. Therefore, by performing the last correction at an internal point of the grid, 
a more accurate Isolution (that is, smaller truncation error for a given order) is obtained. 

The second (and not so intuitive) reason is that the introduction of back corrections induces 
numerical stability. It is well known that the general stability boundaries for the predictor- 
corrector methods of the Stormer-Cowell type decrease geometrically with increasing order. 
That is, when the higher order Stormer-Cowell methods are used on a practical problem, 
the step size is severely constrained by stability considerations. This is not the case with the 
methods using back corrections. In fact, in some cases, these methods possess general stab- 
ility boundaries which are 30 to 40 times larger than those of the Stmmer-Cowell methods of 
the same order. Moreover, the stability regions of these methods do not exhibit geometrical 
decay but remain relatively unchanged up to the eighteenth or nineteenth order. 
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The back<orrection methods are given by: 

k 

■•‘2) ^ -(m + 1) ^ or. f . 

^♦1 ^ ' n«m ' ' n-iii>>l j n-j 

j=0 

and 

k 

^+i-« = (n» + 2 - 8) x„.^ -(m + 1 - 8) Pjb 

j=o 

for 8 = 0, 1 , .... m. The position vectcnr of the satellite at time t^ is denoted by , and 
the acceleiati(Hi vector, by f^ . The numbor of bade coirections is m. I'he coefficients |a^ 
and {0.J are determined in such a way that the highest possible order is obtained for a 
given k. It should be noted that for m = 0, the equations reduce to those of the Stormer- 
Cowell method. Therefore, the back-correction methods are simply a generalization of the 
classical Stormer-CoweU method. 

These methods can be used with pseudo-evaluations in the following algorithm, provided 
the acceleration can be separated according to dominant and perturbing terms: 

a. Predict a value, . for ^ i using equation 1 . 

b. Evaluate the dominant and perturbing accelerations using this value for ^ , saving 
the perturbing acceleration for subsequent calculations. 

c. Set 8 = 0. 

d. Using equation 2, obtain a corrected value, x® * ’ j , for , g . 

e. Reevaluate only the dominant acceleration using x® * ‘ .g and obtain f„ + , .g by 
adding the previously calculated perturbing acceleration at t^^ , g. 

f. If 8 = m, proceed to the next step of the integration. Otherwise, set 8 = 6 + 1 and 
go to (d) above. 

In problems involving earth satellites, the forces can readily be separated according to dom- 
inant and pertiubing terms. In such problems, these algorithms are especially efficient, 
because all evaluations of the forces after the first are simply pseudo-evaluations requiring 
only the reevaluation of the dominant forces (in this case, the two-body force). 

The methods using back corrections have been tested on several problems. A numerical 
int^ation of the orbit of the Applications Technology Satellite-F (ATS-F) has been per- 
formed in which case the errors were smaller by two or three orders of magnitude when 
compared to the classical Stormer-Cowell method using pseudo-evaluations. 

The results shown in table 1 are obtained when the seventeenth-order back -correction 
algorithm is applied to the to the Geodetic Earth Orbiting Sate!Iite-C (GEOS-C) orbit. 


( 1 ) 


( 2 ) 
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Also shown are the results obtained using the twelfth-orda' Stormer-Cowell algorithm 
presently avaflable in the Goddard Trajectory Determination System (GTDS). The errors 
are obtained by comparing with the results obtained using a step size of 40 seconds. It 
should be noted that for reasonable step sizes, the seventeenth-order back-correction method 
provides greater accuracy than the twelfth-order Stormer-Cowell method. Since pseudo- 
evaluations are used in both methods, the computer time required per step is the same 
(within a few percent) fOT both methods. 


Table 1 

Numerical Results for GEOS-C Satellite 


Step Size (s) 

Errors in Position after 24 Hours (km) 

\2th Order 
Stormer-Cowell 

nth Order 

Back-Correcticm Method 

60 

5.7 X 10 ® 

1.9 X lO-'' 

80 

3.2 X \(T* 

8.8 X 10^“' 

100 

1.1 X 10^^ 

1.5 X 10-^ 

200 

6.9 X 10-' 

2.2 X 10-* 


The stability region of the twelfth-order Stormer-Cowell method is, in fact, slightly smaller 
than that of the seventeenth-order back-correction method. A higher order Stormer-Cowell 
method would no doubt have exhibited greater accuracy than the twelfth-order Stormer- 
Cowell method, but the stability region would have been considerably reduced. Consequent- 
ly, the largest meaningful step size attainable with such a higher order method would have 
been much smaller. 

The methods using back corrections appear to be more efficient than the classical methods 
for problems in which the dominant and perturbing forces can be readily separated and in 
which the evaluation of the perturoing terms requires much more computer time than the 
evaluation of the dominant terms. 
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SOLUTIONS OF THE MOTION OF SYNCHRONOUS SATELLITES 
WITH ARBITRARY ECCENTRICITY AND INCLINATION 


Paul E. Nacozy and Roger E. Didil 
University of Texas 
Austin, Texas 


ABSTRACT 

A first-order, semianalytical theory for the long-term motion of resonant satellites is presented. 
The theory is valid for all eccentricities and inclinations and for ail commensurability ratios. 
The method aUows the inclusion of all the zonal and tesseral harmonics as well as luni-solar 
perturbations and radiation pressure. 

The method is applied to a ^nchronous satellite including only the Jjand harmonics. 
Global, long-term solutions for this problem, eccentricity, argument of perigee, and inclina- 
tion are obtained. 

INTRODUCTION 

The method of solution presented here is based on a modified Von Zeipel method applied 
to resonant satellite systems with three degrees of freedom. The method allows any in- 
clination and eccentricity and is a first-order semianalytical theory, yielding only the long- 
term perturbations. 

The procedure first eliminates the short periodic terms numerically by a classical Von Zeipel 
averaging process. Since the averaging is performed numerically, developments in the 
eccentricity and inclination are not necessary. For the same reason, it is straightforward 
to add luni-solar perturbations and radiation pressure. 

After the first Von Zeipel transformation, we have an averaged Hamiltonian that is a tabulated 
function of the variables: argument of perigee, longitude of the node, eccentricity, and in- 
clination. The averaged Hamiltonian is obtained as a numerical table, not m analytical 
function. 

The next procedure is then to eliminate the angular variable corresponding to the resonance 
argument, or the critical argument, t . a modified Von Zeipel tiansformatiort. The principle 
here is that the averaged Hamiltonian at resonance is a minimum with respect to the critical 
argument. This principle was introduced by Hoii in 1960 for his satellite critical inclination 
theory. It was used in a slightly different form by Musen and Bailie in their satellite tesseral 
resonance theory in 1 962. During the transformation, the new Hamiltonian is set equal to 
the value of the old Hamiltonian at its minimum. 
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After these two canonical transformations, a one-degree-of-freedom system remains, and 
in that system the Hamiltonian is a constant. We are now able to analyze the new one- 
degree-of-freedom Hamiltonian and obtain all long-period and secular perturbations, without 
having to develop a %ries in eccentricity or inclination. 

The procedure and algorithm of the modifled Von Zeipel method that we are using was 
developed by G. Giacaglia in 1965 (and described in detail in 1969). He applied the method 
to resonant asteroids in a circular restricted three-body problem in 1968. It was then applied 
to the elliptic restricted problems in 1970 by Giacaglia and Nacozy. The results of the 
elliptic restricted problem are in agreement with the solutions of asteroidal motion performed 
by Schubart (1968) using methods of numerical averaging. 

The method was applied to the Huto-Neptune system by Nacozy and Diehl (1974) and the 
results agree very closely with the numerical integration performed by Williams and Benson 
(1971), who used methods of numerical averaging over a four and one-half million year 
period. 

Recently, we have adapted the method to satellite systems in resonance with the tesseral 
harmonics and have applied it to the -I] *'^22 Problem. The adaptation and application are 
presented below. 

DESCRIPTION OF THE METHOD 

The method uses the Delaunay equations of motion. The angular variables are modified 
slightly so that one variable is the mean anomaly, V^. A second variable is the critical 
argument tor the synchronous satellite, Y, , and a thud i: the argument of perigee, Yj. 

Letting the X variables refer to the conjugate momentums, and where F is the Hamiltonian, 
we have 


Xo = L-H Y„ = K 

X, = H Yj = 6 + w + n+e 

Xj=G-H Yj=w 

F = 1 /2 • L* + H + geopotential (+luni-solar forces) 

The first canonical transformation is the classical Von Zeipel transformation where the new 
Hamiltonian, F*, is defined as the average of the old Hamiltonian over the mean anomaly 
from 0 to 2jt: 


Jir 

F*(Xo*.X,*,X,*;-,Y,*.Y/) - f F(X,, X,.X,; Y^, Y,, Y,)dy^,. 
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After the averaging, we replace the X and Y variables by X* and Y* and then a new Hamil- 
tonian is obtained, free from short periodic terms. This quadrature is performed numerically, 
not analytically. The quadrature is performed for various values of X^ *, Xj •, Y^ *, and Y^ *, 
and the Hamiltonian is obtained as a tabulated function of these variables. 

The new Hamiltonian, F*, has a minimum at Y^ * near 90°. This minimum corresponds to 
the libration center for Y, (also 270°). The exact value of Yj ♦ depends on the value of Y^ ♦ 
(the perigee). The partial derivative of the Hamiltonian with respect to Yj *, at Yj * near 90°, 
is zero. 

The next canonical transformation defines a new Hamiltonian, F**, to be the value of F* at 
its minimum value: 

F** (X„**, X, *•, Xj**; -, -, Yj**) 

= F* 

Y,*3s90°(or270°). 

This canonical transformation is that used in 1960 by Hori for his critical inclination solution, 
and we are using the same principle here. Since F** does not depend on Y, **, we have: 


X,- = 


3F** 

=0. 

3Y •• 


This gives a quasi-integral for Xj •* : 

Xj •* = H** = constant = >/l~c**2 cosi"*. 


This relation between e** and i** is the component of the angular momentum in the Z 
direction. Since F** is a constant of the motion (conservative system), we now have two 
relations between three variables-eccentricity, inclination, and perigee— and we can rewrite 
the two relations as: 

F**(e,i,<»)) = constant 
Ii**(e, i) = constant 

These two functions deflne curves in an e versus w plane with H** as a parameter or curves 
in an i versus cj plane, also with H** as a parameter. This then gives us the long-term solution 
for e, i, and to. 


DETAILS OF AN APPLICATION 

The first part of the numerical results for the Jj - Jjj problem are shown in figure 1 for F* 
versus the critical argument, Yj *. Only one plot is ^own here for Yj * (or cj*) 90°, but many 
other plots were obtained for different values of Yj * (<o*). We have values of eccentricity 
ranging from 0.0 to 0.9, but only six values are shown in figure 1 . Also, only one value for 
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the parameter H** is given, corresponding to H** to 0.3. It should be noted that there is 
a minimum of F* at Y, * = 90°. This corresponds to the libration center for the synchronous 
satellite at Yj • = 90°. 

We define the new Hamiltonian to have the minimum value of the old Hamiltonian, and then 
we plot the new HamUtonian. This is shown in figure 2, for F** versus Y, **, where Yj •* goes 
from 0“ to 180°. The problem is periodic inir, so this will repeat from 1 80° to 360°. 

The plot of figure 2 may be considered as a three-dimensional plot by visualizing an eccentric- 
ity axis as perpendicular to and coming out of the plane of the paper. This will then produce 
a two-dimensional surface and it can be seen that there will be a valley in the surface centered 
atYj** = 90° and e**= 0.44. 

Since F** is a constant for a given trajectOTy, we can construct planes parallel to the e versus 
cj plane (in the three-dimensional plot), one plane for each value of the constant Hamiltonian. 
In other words, we take several contour levels and each of the levels corresponds to a constant 
Hamiltonian and hence to a certain trajectory. The contour curves give the long-term global 
solutions for eccentricity versus perigee and are shown in figure 3. Cbntour levels near the 
bottom of the valley produce libration. At the bottom of the valley is the libration center, 
and at the top of the valley there is circulation. 

Consider an ec jntricity of 0.1. We see in figure 3 that, for this eccentricity, there is cir- 
culation of the perigee by following the curve, beginning at e = 0. 1 and Y^ ** = 0. In other 
words, the e** versus o)** on that curve corresponds to an actual trajectory, the long-term 
solution to the problem. 

If we have an eccentricity of about 0.4, we have libration of the nerigee about the value of 
perigee equal to 90°. At nearly e = 0.44, there is a stationary solution, and this corresponds 
to a periodic solution. Apparently, these periodic solutions have not been found prior to 
this work for the Jj - Jjj problem. 

Using the quasi-integral relation, we can obtain corresponding plots for the inclination versus 
the perigee, and this is shown in figure 4. The libration center that was shown in figure 3 
corresponding to c = 0.44 shows up at about i** = V'0° and <o** = 90°. 

The stationary solution (libration center) that we have found and presented here has an 
eccentricity of 0.44, an inclination of 70°, and a perigee of 90°. We have found that as H** 
is varied, these stationary solutions form a family of periodic solutions. 

CONCLUSION 

Our future studies will be to trace the family of periodic (stationary) solutions. We will also 
add more zonals and tesseral harmonics, the effects of the sun and the moon, and radiation 
pressure to determine what effect they will have on the family of stationary, librating, and 
circulating solutions. 


178 




1 



180 



REFERENCES 

Giacaglia, G.E.O., 1965, Doctoral Dissertation, Yale University. 

Giacaglia, G.E.O., 1968, Smithsonian Astronomical Observatory Special Report No. 278. 
Giacaglia, G.E.O., 1969, Astron. J., 74, p. 1254. 

Giacaglia, G.E.O. and P.E. Nacozy, 1970, appears in Per: .die Orbits, Stability and 
Resonances, ed. by G. Giacaglia, D. Reidel Publishing Co., p. 96. 

Hori, G., I960, Astron. J., 65, p. 291. 

Musen D. and A. Bailie, 1962,/. Geophys. Res., 67, p. 1 123. 

Nacozy, P. and R. Diehi, 1974, Bull. Astron. J. , 6, (2). 

Schubart, ?., 1968,/lrtron, /., 73, p. 99. 

Williams, J. and G. Benson, 1971, Astron. /. , 76, p. 1 67. 

DISCUSSION 

VOICE: What do you mean by stable stationary solutions? 

NACOZY: I mean stability in the sense that if you place a particle on any one of the contour 
curves (corresponding to a trajectory), any slight deviation from the curve will cause the 
particle merely to move to another adjacent curve. 

VOICE: What do you mean by a <«Masi-integral? 

NACOZY: The integrals that I h ; presented are accurate only to first order and hence 
are often re"'rred to as quasi-integrals. 
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SYSTEM DESIGN IMPACT OF GUIDANCE AND 
NAVIGATION ANALYSIS FOR A SEPS 1979 
ENCKE FLYBY 


Philip Hoqg 

Martin Marietta Corporation 
Denver, Colorado 


This is a report on a study that Martin Marietta did for Rockwell, who in turn was doing a 
feasibility study for the solar electric propulsion stage (SEPS). The emphasis was on system 
feasibility, and we tried to merge the guidance and navigation requirements into the total 
system. The primary emphasis was the 1979 Comet Encke flyby, which at that time was 
the first proposed SET'S mission. We also looked at some system requirements fw other 
missions. 

Many of the system parameters are affected by guidance and navigation requirements. The 
most important ones are thrust control authority (that is, how much additional control is 
needed in the thrust subsystem to implement trajectory ctnrections), thrust performance 
tolerances, thrust vector control, propulsion time, the additional time required for adjust- 
ing the trajectories, fuel requirements, guidance updates, and the types of earth-based 
navigation and communications needed. Other parameters affected are the onboard naviga- 
tion subsystem (how good must it be, if it is indeed needed) and the type of trajectory and 
terminal errors that occur (that is, control and knowledge). Control is the dispersion of the 
actual about the reference, and knowledge is the dispersion of the estimated about the 
actual. 

The baseline mission for this particular stage was a launch in March 1979 and encounter in 
November 1980. Figure 1 shows an ecliptic projection of the flyby. 


L - LAUNCH POSITION 
C - tNCOUNTCR IL * S«2 DAYS) 



Figure 1. Ecliptic projection of the 1979 earth-to*Encke flyby. 
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One of the first things we did was to take the mission analysis results that Rodcwell gen- 
erated and produce a more realistic, optimized and targeted trajectory in which we imposed 
reasonable control policies. Listed below are the mission data: 

• Launch date-March 25, 1979 

• Airivaldate— November 7, 1980 

• Launch VHE— 7. 1 8 km/s 

• Initial mass-1988 kg 

• Initial power— 21 kW 

• Housekeeping power-0.650 kW 

• Thruster effic’ency— 64 percent 

• Propulsion time— 523 days 

• Coast time— Initial, 64 days; final, 5 days 

• Arrival VHP— 3.16 km/s 

• Arrival RCA— 1 100 km 

• Arrival mass-1456 kg 

We employed constant cone, clock, and thrust over fixed time segments, and we imposed 
whatever constraints were necessary; for example, the final coast time of 5 days was the 
minimum acceptable for science. There was also an initial coast and a total thrust time of 
523 days. We arrived at 3 km/s at 1100 km closest approach. The encounter was 30 days 
prior to perihelian for this particular study. 

The baseline guidance and navigation strategy-and we are talking mostly about the approach 
phase now— assumed simultaneous or continuous coverage from three Deep Space Network 
(DSN) stations over tte last 40 days. We assumed an onboard optical system ^ilar to 
Mariner-10. Optical or onboard observations were taken twice per day starting at 30 days 
prim to. encounter, and each optical measurement contained Encke and three identifiable 
stars. We estimated the vehicle state, thrust biases, and the Encke ephemeiis, and we con- 
sidered the process noise that is generated by the thrusters. Guidance u\ 'ates were per- 
formed every 4 days, and we assumed that we could control the trajectory just by biasing 
the nominal thrust controls. The baseline strategy was not intended to meet all the system 
requirements, but it was our first try at it. 

The dynamic and measurement error sources that we assumed are listed below: 

• Dynamic error sources 

a. Launch errm — Position, 3 km; velocity, 5 m/s; mass, 1 kg 

b. Thrust bias — Magnitude, 2.2 percent; direction, 0.035 rad 



c. Thrust noise — Magnitude, 3.S percent; correlation time, S days; direction: 
0.010 rad; correlation time, 3 hours 

d. Encke ~ Position, 10,000 km; velocity, lOOOkm/day 
• Measurement error sources 

a. DSN station location — Spin radius, l.S m; longitude, 3.0 m; height, 1 0.0 m 

b. Dopplo- noise — Two-way, 1 mm/s; three-way, 0. 1 mm/s 

c. Range noise —Two-way, 3 m; three-way, 10 m 

d. Optical — Resolution, 30 arc-seconds; center flnding, 10 km 

SEPS is launched with the Titan^Tentaur, so we had typical insertion uncertainties. The 
thrust uncoiainties consisted of both bias and noise. We estimated the bias and the noise, 
whidi has a time-varying component, hence the correlation times associated with them were 
considered. The values listed are for a sin^e thruster; there are eight operating thrusters for 
this vehicle. The a priori error for Encke was assumed to be very pessimistic, although at 
that time it was a reasonable uncertainty of 10,000 km and 1000 km per day, or about 10 
m/s. 

TyiHcal station location uncertainties were assumed. Since we used simultanemis range and 
range rate data from earth to the spacecraft, we also induded three-way uncertainties. Our 
onboard optical uncertainties consisted of an error or 30 arc-seconds nmse and an uncer- 
tainty of 10 km to approximate the uncertainty between the comet’s center of brightr 
and the center of gravity. 

One of the first things we studied was the approach geometry, to see if we could get an idea 
of what was happerung. Figure 2 is a view of the last 30 days prior to encounto-. It can be 
seen that the thrust vector is almost retrograde, which means t^tat it is on a flat trajectory 
relative to Encke, and there is very little curvature. There is also the possibility of the 
thrust plume into'fering with any instruments that are sensing as Encke is approached. 
However, cutoff occurs at 5 days prior to encounter. The position of the earth is changing 
very rapidly both in direction and in declination, which would indicate that earth-based 
tracking of the spaceoaft would be very good. On the other hand, because of the flatness 
of the trajectory, the onboard optics probably will not be as good. 

Figure 3 is a plot of the inertial position uncertainties due to the estimation of both Encke 
and SEPS. Because simultaneous data are employed, we ^t near-ballistic results for the 
spacecraft. However, our onboard optics just barely get below the a pnori level. 

Most of this is along-track error. However, a 6000-km along-tradc error is about 1 0 minutes 
uncertainty, and that would seriously affect pointing and slewing rate* Because the stage 
and Encke are dynamically uncoupled, the Encke relative uncertainty is basically the RSS 
of these two, which is then dominated by the Encke uncertainties. 


184 



7 

S 

$ 

3 
2 
1 
0 

012 345 0709 10 

no*KM» 

Figuie 2. Encke ^lyby relative approach geometry for last X days prior to encounter. 

As shown in Iflgure 4, the velocity shows only fair improvement. It is stiU not very good, 
either that of Encke or of the stage. 

Snce we found that the Encke relative uncertainties are dominated by the a priori ephemeris 
error, we looked at a number of different a prioris and found that Encke relative uncertain- 
ties are a priori sensitive, pr’marily in the velocity component (figure 5). For example, we 
found that, if we have velocity uncertainty, we will get very good estimation errors. 
However, if we ha^-e any reasonable amount of a priori velocity uncertainty, the optics 
cannot comper U, which means that the comet’s relative uncertainty remains unchanged. 

We did takf some covariances from Bob Farquhar and used them as our a priori. They looked 
more lik' >ur zero velocity case, although not quite as good, or about 400 km terminal 
uncer*.;inty. 

F . ae 6 illustrates the magnitude of the thrust guidance corrections, that is, biases to the 
nominal thrust control policy. The biases never exceed more than 3 percent in thrust mag- 
nitude and are ir-egularly shaped because of our strategy. The curve could have been smoothed 
out by pTrin . the guidance updates at different intervals and by employing a different 
control p )iicy. But in magnitude, the bias is not more than 3 percent, and in pointing, it 
islc.< . man 2°, lo. 
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Figure 5. Encke ephemeiis uncertainty. 

Figure 7 illustrates the terminal control error; that is, the error obtained after all the orbit 
determination and guidance have been employed. The shaded area is the inertial spacecraft 
uncertainty, and the total bar is the Encke relative uncertainty. We chose 1000 km as a 
guidance success zone, which is our maximum acceptable error for a successful mission. 
Obviously our baseline missed that by quite a bit. If we look at the effect of a priori 
ephemeris uncertainties, we could get within the success zone if we have zero velocity un- 
certainty in Encke. However, that is not a realistic case. If we assume some currently reason- 
able error, we will probably just barely make the guidance success zone. Because of the 
dominance of the ephemeris uncertainty, even if we eliminated earth-based ranging of the 
stage, the Encke relative uncertainty would not be affected, but the stage uncertainty would 
be. If we ever do reduce the ephemeris error down to some low level, we wQl still need 
earth-based ranging, and the same thing goes for simultaneous range and range rate data 
from the earth. 

The conclusions for this particular study are as follows: Cruise guidance and navigation 
requirements (not discussed) were minimal compared to the approach phase. The Encke 
relative approach error was dominated by the ephemeris uncertainties, particularly the 


187 




MAGNITUDE (XI 



Figure 6. Encke flyby thrust guidance corrections. 
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velocity components. The earth-basv,ci tracking was very good; however, quality three-way 
data are still needed from the DSN b'utions. Optical navigation was only fair, due primarily 
to the flat trajectory. Finally, the thrust control authority, that is, the thrust updates, were 
within acceptable tolerances. 

Table 1 lists some other missions which might drive the design of this stage, which is supposed 
to be a multimission spacecraft, and we identified a few missions which might drive the 
g' lidance and navigation subsystem design. One was the Encke flyby, because that was the first 
proposed SEPS mission. The next was an earth-orbital mission, because of its unique char- 
acteristics. The Encke rendezvous was chosen primarily because it might place great demands 
on an onboard sensor and also because it was the second interplanetary mission. The 
Mercury orbiter was chosen because it had a high thrust acceleration at approach, which 
meant that the process noise in the thrust would dominate this particular mission. Finally, 
we chose an outer planet mission, or an outbound mission, in this case a Phobos/Deimos 
rendezvous. The data in the table are just guesses as to what the impact might be on these 
subsystems. 


Table 1 

Mission/Subsystem Impact* 


Mission 

Launch 

Flight 

Time 

(days> 

Approach 

Guidance 

Thrust 

Thrust 
Vector and 
Attitude 
Control 

Earth- 

based 

Tracking 

Operations 

and 

Communications 

1 

Science 

Encke Flyby 

1979 

630 

M 

L 

t 

M 

L 

H 

Earth Orbital 

1980 

50-1001 

H 

H 

H 

H 

H 

- 

Encke Rendezvous 

1981 

1100 

H 

H 

M 

M 

M 

H 

Mercury Orbiter 

1984 

400 


H 

H 

H 

M 

M 

Phobos/Deimos 

Rendezvous 

1984 

260 

H 

1 

M 

I 

M 

M 

L : 

H 

H 


*(L = low, M = moderate, H = high). 
tPcr one-way trip. 


If we dismiss ephemeris uncertainties for the time being, the limiting factor is the thrust 
uncertainty; that is, the noise or uncertainty in the performance of the engines. Figure 8 
shows the closest approach uncertainty. The Encke flyby is relatively insensitive to thrust 
magnitude uncertainties, primarily because of the good tracking from the earth, which can 
minimize the effect of the thrust uncertainties. The Encke rendezvous in 1981, which is a 
1984 encounter, does not do as well, because the earth geometry relative to the spacecraft 
and Encke is not as good. The Mercury orbiter, because of the high acceleration at approach, 
is ;'ery sensitive to thrust error. Of course, this means that if we have a 1000-km guidance 
success zone, we had better reduce the thrust errors considerably for the Mercury orbiter. 
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Figure 8. Sensitivity of encounter error to thrust error. 

The following recommendations were made and will probably be conducted in a foPow-on 
study: 

• Improve the small body, both comet and astroid, ephemeris determination. Part 
of this involves integrating the earth-based telescope observations of Encke with 
the DSN and onboard optics measurements. 

• Reduce thrust noise level, either through hardware changes or through better 
orbit determination and measurement. 

• Continue development of the simultaneous and/or differenced data types, that is, 
the quality earth-based three-way data. 

• Investigate the impact of other missions and combine their requirements into a 
common spacecraft. 

• Study alternate mission strat^ies. For example, since we do have thrust at the end, 
we might shape the trajectory and possibly get some curvature to minimize the 
along-track error. 
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GROUND TRUTH APPLICATIONS TO ORBIT REFINEMENTS 


Robert L. White 

The Charles Stark Draper Laboratory 
Cambridge. Massachusetts 


During the past few years, the C. S. Draper Laboratory has been performing various studies 
for Goddard Space Flight Center pertaining to the determination and use of spacecraft 
attitude and orbital ephemeris data to improve the mapping accuracy of an earth-observing 
multispectral scanner. These studies were conducted for an Earth Observation Satellite 
(EOS) that is assumed to be in a circular, sun-synchronous orbit with an altitude of 1 000 
km. 

At present, an investigation is being made into the use of known ground targets (that is, 
landmarks) in the earth sensor imagery, and also stars in combination with known ground 
targets, to estimate the spacecraft attitude, orbital ephemeris, and the bias drifts of three 
strapdown gyros. The present study is a covariance analysis where both the Kalman filter 
and Fraser two-filter smoother are used to process star and landmark measurements to obtain 
a statistical indication of performance. Star measurements are used to update attitude and 
gyro bias drift (that is, 6 parameters), and landmark measurements are used to update all 
12 state parameters. This study is, for the time, restricted to the use of landmarks in the 
continental United States, Alaska, and Hawaii, since the primary interest in spacecraft 
attitude and orbital ephemeris is assumed to be during the observation passes over these 
regions. 

The geometry of a typical pass over the United States is shown in figure 1 . The spacecraft 
maintains a local vertical orientation as shown in the figure, where the body axis, Zg , is 
always directed toward the subsatellite point. The star tracker is assumed to be a body-fixed 
instrument whose optical axis is directed toward zenith. As the spacecraft circles the earth, 
the stars pass through the 8° square field-of-view (FOV) and are electronically tracked. For 
the purposes of this study, these stars are artificially generated with random positions in the 
FOV at the times of measurement. These measurements are uniformly distributed through- 
out the orbit and only one measurement is made on each star as it passes through the FOV. 

On board the spacecraft, there is assumed to be a ...ultispectral scanner whose beam is 
directed downward and sweeps back and forth across the ground track to generate a swath 
of imagery 145 km (90 miles) wide (see figure 1). In the present study, a landmark meas- 
urement represents the line-of-sight (LOS) of the scan beam at the time of landmark 
observation. This LOS is completely defined in body coordinates by the scan beam angle 
for which a random value (within ±4.8°) is selected for each landmark measurement. 
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Table 1 defines the five different landmark observation cases used in this study. The case 
numbers denote the number of orbital passes made over the regions of interest. A pictorial 
representation of case 4 is shown in figure 2, where the plus signs (+) indicate the nominal 
locations of the landmarks used for update purposes. 

The nominal conditions used to generate mo.,i of the performance results are listed below; 

• Orbit - Sun-synchronous circular orbit with an altitude of 1000 km 

• Gravity model — Central force field 

• Spacecraft attitude history — Local vertical with rotation only in pitch 

• Star measuicments — Error: 5 s ( lo)/axis; number: 20 per orbit (evenly distributed) 

• Lanumark measurements — Position. error: 15 m ( 1 a) in downrange and crosstrack; 
number: two/pass 

• Gyro enor - Random drift rate (lo) = 0.01°/hr (white noise; quantization (la) = 

0.1 s 

It can be seen that relatively simple models were adopted for the spacecraft attitude history 
and the gravitational model, since it was felt that these would be sufficient for the purposes 
of this investigation. Previous studies of attitude determination have shown that the perform- 
ance results for a nominal attitude history are in fairly close agreement with those for an 
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Figure 2. Landmark observi».ton m$<- ' 
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Tabic 1 

Landmark Observation Cases 


Case 

Number 

Pnss 

Number 

North Latitude (deg) 

Region* 

At Start 

At End 

1 

1 

50 

30 

USA 

2 

1 

50 

30 

USA 


2 

50 

30 

USA 

3 

1 

50 

45 

USA 


2 

50 

30 

USA 


3 

50 

45 

USA 

4 

1 

50 

30 

USA 


2 

50 

30 

USA 


3 

65 

60 

Alaska 


4 

65 

60 

Alaska 

5 

1 

50 

45 

USA 


2 

50 

30 

USA 


3 

50 1 

45 

USA 


4 

65 

60 

Alaska 



and at 

20 

Hawaii 


5 

65 

60 

Alaska 


*USA denotes oontinental USA. 


attitude history that deviates from nominal by the amount anticipated in EOS. Deviations 
in attitude due to the various disturbing torques can be accurately indicated by the space- 
(7aft gyros. A central force Held was used for the gravity model. It was felt that the inclu- 
sion of J, and the higher gravitational harmo*- v vould add undue complexity to the prob- 
lem without shedding any additional light on thw merits of using star and landmark measure- 
ments. With the possible exception of Jj , the inclusion of ^anr- itiics does not 

produce a significant change in the geometry of the ba.«’' ' ' ■; 'he period of 
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interest (a few orbits). Consequently, if the higher harmonics were included in the simula- 
tion and were accurately accounted for in the {M'opagation of ephemeris data (that is, no 
harmonic uncertainties), it is felt that the performance wouM be somewhat the same. With 
regard to the existence of uncertainties in these harmonics, it is felt that this is a problem 
which all techniques must face. 

The nominal landmark measurement error adopted for this study was IS meters (lo) in 
downrange and crosstrack. This value does not represent the most recent estimate of what is 
anticipated for the proposed EOS; it merely represents the anticipated uncertainty in 
establishing the location of a ground control point (landmark) on the earth by other means 
(for example, surveying). In a mote realistic situation, we wouki also include the error in 
determining the scan beam angle at the time of observation and also the errors associated 
with image resolution and the method used to identify landmarks in the imagery. Since 
most of these other sources of error had not been firmly established for B3S at the begin- 
ning of this study, they were not considered when adopting the present nominal value. 

This was felt to be an acceptable approach, since the plan was to generate sufficient sensitiv- 
ity data to show the effect of using different values of the impo'tant error sources and 
parameters. 

The nominal values used for the initial state uncertainties are as follows: 

• Attitude (pitch, roll, yaw) — 60 arc-s (each) 

• G)to bias drift — 0.03“ /hr (each) 

• Ephemeris position 

a. Altitude — 20 m 

b. Downrange — 50 m 

c. Crosstrack — 20 m 

• Ephemeris velocity 

a. Altitude — 0.05 m/s 

b. Downrange — 0.02 m/s 

c. Crosstrack — 0.02 m/s 

In figures 3 and 4, the performance in estimating spacecraft position is shown for both the 
Kalman filter and the Fraser two-fUter smoother. The results are for landmark observation 
case 4. This case, like all others, was initiated at the ascending node of the orbit and 
completed one pass over the north polar region before makmg the first pass over the United 
States. It is seen that the Kalman filter does an effective job in reducing the position un- 
certainties after two passes over the continental United States, and with only two landmark 
updates per pass. The significant improvement in filter performance after two passes is 
primarily due to the large reducticm in the velocity uncertainties on the second pass. Both 
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Figure 3. Kafinan fitter and Fraser smoother downrange estimation uncertainties 
for landmark observation case 4. 




Figure 4. Kalman f iUer and Fraser smoother estimation uncertainties for crosstrack 
and attitude for landmark observation case 4. 
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figures 3 and 4 show that the smoother does a much better j<^ than the filter for those por- 
tions of the orbit away from the United States. It should be noted that each smoother curve 
exhibits somewhat the same minimum during each pass over a landmark update region, 
while those for the filter are reduced from pass to pass. It is also noted that the filter per- 
formance eventually approaches that of the smoother, and it is for this reason that the fflier 
was used in place of the more complicated smoother to generate most of the performance 
res'ilts of this study. 

In table 2, the performance in estimating aU 1 2 state parameters is shown for both the filter 
and the smoother. The filter data represents the uncertainties at the end of the last pass, 
while that fm* the smoother is for a time at the end of the first or second pass. Note that 
the attitude performance is very good This is primarily due to star updates. The reason 
that the yaw performance is not as good as that for pitch and roll is due to the fact that very 
little yaw information is obtained directly from the star tracker because die stars are dose 
to zenith. 

Table 3 illustrates the performance when either landmark or star updates are not used. The 
first set of data is a repeat of the nominal performance of table 2 (using star and landmaric 
updates) and is shown here for purposes of comparison with the other two data sets. The 
second set of data represents the performance when only landmark updates are used. It is 
seen that the uncertainties in pitch and dowtuange position (designated “range” in the table) 
continue to grow as one goes from landmark observation case 1 to case 5. This dearly 
indicates a lack (rf* observabUity in this approach. It is also interesting to note that the down- 
range uncertainty in meters is numerically about five times larger than the pitch uncertainty 
in arc-seconds. Since 1 arc-s subtends about 5 m at a distance of 1 000 km, the results indicate 
that the Alter, after overcoming the initial state uncertainties, ends up applying the same 
equivalent update to both pitch and downrange position. The landmark measurements 
[H'ovide very little information to distinguish between the two state parameters. Consequent- 
ly, the existing uncertainties in gyro bias drift and spacecraft velocity cause the fatch and 
dowiuange position uncertainties to grow with time. It should be noted that the same 
numerical relationship (5 to I) occurs between roll and crosstrack; however, the uncertain- 
ties in these parameters do not grow with time since they are bounded by nature. On the 
basis of the pitch and downrange perforiiiance, it would theref^/re seem that there is no 
useful purpose to be gained in using only landmarks. However, it has been found that very 
strong nigative correlations do occur between the errors in pitch and downrange position 
and also between the erron in roll and crosstrack position. These correlations are such as to 
greatly nullify the effects of these errors on a mapping process that makes use of attitude 
and ephemeris data updated with only landmark measurements. 

The third set of data in table 3 shows the performance when only star updates are used. 

These data were generated as a matter of interest, since star measurements can only be used to 
update attitude and gyro bias drift. It is seen that the attitude performance is almost as 
good as that of the first set. The uncertainties shown for the spacecraft position components 
represent the natural growth of these quantities, and it is seen that the downrange position 
uncertainty grows more rapidly than that in the second data set. 
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Table 2 

Kalman Filter and Fraser Smoother Estimation Uncertainties 
for Different Landmark Observation Cases 




State Estimation Uncertainties (la) 

Landmark 

Observation 

Case 

Filter (F) 
or 

Smoother 

Attitude 

(arc-s) 

1 

Gyro Bias Drift 
(lO-Vhr) 

Position 

(m) 

Velocity 

(m/s) 


(S) 

Pitch 

Roll 

Yaw 

X 

Y 

Z 

Alt, 

Range 

Track 

Alt. 

Range 

Track 

Initial State Uncertainties 













(AURuns): 


60 

60 

60 

30 

30 

30 

20 

50 

20 

0.05 

0.02 

0.02 

1 

F 

3.1 

3.1 

20 

20 

2.1 

4.0 

41 

21 

14 

0.05 

0.03 

0.02 


S 

2.5 

2.3 

20 

20 

2.1 

4.0 

30 

18 

14 

0,04 

0.02 

0.02 

2 

F 

1.9 

1.6 

14 

14 

0.4 

1.0 

14 

15 

10 

0.02 

0.01 

0.02 


S 

1.2 

1.5 

14 

14 

0.4 

1.0 

14 

13 

10 

0.01 

0.01 

0.02 

3 

F 

1.5 

1.4 

12 

12 

0.2 

0.7 

14 

12 

8 

0.u2 

0.01 

0.02 


S 

0.8 

1.3 

12 

12 

0.2 

0.7 

14 

8 

8 

0.01 

0.0 1 

0.02 

4 

F 

1.3 

1.3 

9.8 

9.8 

0.1 

0.6 

12 

11 

7 

0.02 

0.01 

0,02 


S 

0.8 

1.1 

9.7 

9.8 

0.1 

0.6 

12 

7 

7 

0.01 

0.01 

0.02 

5 

F 

1.2 

1.3 

8“ 

8.7 

0.1 

0.5 

9 

10 

7 

0.01 

0.01 

0.01 


S 

0.8 

1.1 

8.7 

8.7 

0.1 

0.5 

9 

7 

6 

0.01 

0.01 

0.01 


Note: Smoothed data are for the time at the end of the second pass over the United States except for cases 1 and 2 where 
data are for the time at the end of the first pass. Filter data are for the time at the end of the last pass. 




Table 3 

Kalman Filter Performance for Cases With and Without Landmark or Star Updates 


Landmark 

Obseivation 

Case 

STATE ESTIMATION UNCERTAINTIES (la)* 

Attitude 

(arc-s) 

Position 

(m) 

Pitch 

Roll 

Yaw 

Alt. 

Range 

Track 

Initial Uncertainties Are: 

60 

60 

60 

20 

50 

20 

Nominal Performance with 2 Landmark Updates per Pass and 20 Star Updates per Orbit 

1 

3 

3 

20 

41 

21 

14 

2 

2 

2 

14 

14 

15 

10 

3 

2 

1 

12 

14 

12 

8 

4 

1 

1 

10 

12 

11 

7 

5 

1 

1 

9 

9 

10 

7 

Performance with 2 Landmark Updates per Pass and No Star Updates 

I 

44 

5 

32 

66 

213 


2 

90 

5 

2"» 

54 

441 

20 

3 

161 

4 

22 

50 

782 

20 

4 

229 

4 

19 

38 

III! 

20 

5 

277 

4 

17 

29 

1342 

20 

1 Performance with No Landmark Updates and 20 Star Updates per Orbit 

1 

3 

4 

24 

101 

291 

20 

2 

2 

2 

15 

101 

747 

20 

3 

2 

2 

12 

101 

1226 

20 

4 

1 

1 

11 

99 

1712 

20 

5 

1 

1 

10 

99 

2243 

20 


*Gyro bias drift and spaceaaft velocity uncertainties are not shown. 
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Table 4 shows the sensitivity of performance to variation in the number of iandmarK updates 
per pass and the number of star updates per orbit fcM* landmark obs-a^tion case 2. It is seen 
that fairly good performance is obtained even with one landmark update po* pass and that 
there is no significant improvement when going frcia two (nominal) to five landmark updates 
per pass. It is also seen in the second set of data of table 4 that some variation can be allowed 
in the number of star updates per orbit without seriously affecting the results. 

Table S gives the sensitivity of performance to variation in the star and landmark measurement 
errors for landmark observation case 2. These data ^ve some indication of the measurement 
accuracies needed in order to obtam a desired level of performance. More recent results 
(not shown) have been generated to show the effect of larger variations in these measure- 
ment errors. 
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Table 4 

Sensitivity to Number of Landmark Updates Per Pass and Number of Star Updates 
Per Orbit for Landmark Observation Case 2 


Updates Per 
Pass or 


State Estimation Uncertainties (la)* 


\ttitude (arc-s) 


Position (m) 

Orbit 

Pitch 

Roll 

Yaw 

Alt. 

Range 

Track 

Initial Uncertainties Are: 

60 

60 

60 

20 

50 

20 

Landmarks** 







1 

1.9 

1.6 

15 

16 

18 

11 

2 (Nom.) 

1.9 

1.6 

14 

14 

15 

10 

3 

1.9 

1.6 

14 

13 

13 

9 

5 

1.9 

1.6 

12 

13 

12 

9 

Stars** 







5 

3.3 

2.5 

19 

13.4 

21 

13 

10 

2.5 

2.1 

18 

13.4 

17 

11 

15 

2.1 

1.8 

16 

13.4 

16 

10 

20 (Nom.) 

1.9 

1.6 

14 

13.4 

15 

10 

25 

1.7 

1.5 

12 

13.4 

14 

10 


*Gyro bias drift and spacecraft velocity uncertainties are not shown, 
**Other type measurement (star or landmark) was nominal. 




Table 5 

Sensitivity to Star and Landmark Measurement Error for Landmark Observation Case 2 


Measurement 
Error ( 1 o) 
(Per Axis) 

State Estimation Uncertainties (lo)* 

Attitude (arc-s) 

Position (m) >> 

Pitch 

Roll 

Yaw 

Alt. 

Range 

Track 

Initial Uncertainties Are: 

60 

60 

60 

20 

50 

,fO 

Star Measurement Error 







2s 

0.8 

0.8 

7 

13.3 

12 

8 

5 s 

1.9 

1.6 

14 

13.5 

15 

10 

10? 

3.6 

2.7 

18 

13.8 

21 

13 

ISs 

5.4 

3.3 

20 

14.1 

29 

16 

Landmark Position Error (Downrange and Crosstrack) 





7.5 m 

1.9 

1.6 

12 

10 

11 

8 

IS m 

1.9 

1.6 

14 

14 

15 

10 

30 m 

1.9 

1 7 

15 

16 

24 

13 

45 m 

1.9 

1.7 

15 

16 

34 

16 


•Gyro bias drift and spacecraft velocity uncertainties arc not shown. 
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ATTITUDE DETERMINATION USING 
DIGITAL EARTH PICTURES 


Lawrence P. Gunshol 
Computer Sciences Corporation 
Silver Spring, Maryland 


We at Computer Sciences Corporation have developed a computer program called PICATT, 
which stands for picture attitude determination. This paper describes the particular satellite 
to which this technique of attitude determination has been applied, describes the method 
of solution, and discusses the results that have been attained using the PICATT program. 

The satellite to which we have applied the PICATT technique is the Applications Technology 
Satellite-3 (ATS-3), which has been operational since 1967. The mission characteristics of 
ATS-3 are as follows: It is in a geosynchronous, circular orbit, with very low eccentricity; 
currently the inclination is on the order of S°. It is spin stabilized at approximately 100 
rpm. The dominant torque on the spin axis is the solar radiation torque, which varies from 
10 to 60 arc-seconds per day, depending upon the time of year. 

Observable misalignments have developed between the ATS-3 principal axis of spin and the 
geometric axis of spin. This phenomenon was originally detected by Westinghouse Corpor- 
ation, udng a graphical technique for processing pictures. When Computer Sciences developed 
the PICATT program for ATS-3, we also detected the misalignments between the principal 
axes and the geometric axes. Due to the fact that there are misalignments between the 
principal and the geometric axes, there is a cross-coupling effect on the spin axis whenever 
an east-west stationkeeping maneuver is performed. ATS-3 is positioned at 70° west longi- 
tude, and east-west stationkeeping maneuvers are performed at approximately 3-month 
intervals; whenever one is performed, the attitude is perturbed. 

The Soumi camera has a variable elevation that covers an approximately 18° field of view. 

At synchronous altitude, the earth subtends an angle of approximately 17.4°, and ATS-3 
is able to v»ew the entire earth when it is oriented properly. The elevation angle range covers 
9° above the plane normal to the spin axis to 9° below the plane normal to the spin axis. 

The resolution of the camera, when it is pointing directly at the local vertical, or nadir, is 
approximately 3.7 km (2 n.m.). The spectral bandpass is 4720 A to 6300 A; hence, it is 
a visible system. The analog video signals are processed to form an image on a film, and the 
National Oceanic and Atmospheric Administration (NOAA) produces prints of the earth 
image. ATS-3 originally had the capability of taking a color photograph and was the first 
satellite to take a color photograph of the earth from synchronous altitude. Currently, the 
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capability exists only for taking black-and-white, spin-scan, cloud-cover pictures. It is also 
possible to record the video signals onto a digital tape; these processed digital tapes are used 
for our attitude determination efforts. 

The attitude control requirement for ATS-3 is to maintain the spin axis attitude to within 
1° of orbit normal. More precisely, NOAA desires to take pictures of the earth 2.5 hours 
prior to high noon and 3.5 hours after high noon, for a total elapsed time of 6 hours, or 
one-quarter of an orbit. NOAA also wants to see 55° north latitude and above during this 
time interval. 

As shown in flgure 1 , we use the standard celestial inertial coordinate frame, where the Xg 
axis points toward the vernal equinox and points toward the north celestial pole. Right 
ascension is measured in the conventional sense, counterclockwise from Xj^ . Declination 
is measured positive above the equatoWal plane. 


Ze 



(•) 


^ (SPIN AXIS) 
i 



! 


f 

(bl 


Figure 1. Coordinate systems, (a) Spin axis attitude angles; (b) prindnal and geometric axes. 


As mentioned before, there are misalignments between the ATS-3 princip.’! and geometric 
axes. These misalignments are described by two rotation angles: A geomt trie system is 
defined whereby the X geometric axis points along the symmetric axis of the spacecraft; 
the Z geometric axis points toward the plane of the camera on board ATS-3. We define the 


204 



offset between the principal and geometric axes by first rotating around the Z principal axis 
through the engle ^ (which we call a skew misalignment angle) and then rotat- ^ about the 
Y geometric axis through the angle d (which we describe as an offset misalignment angle). 

As seen in figure 2, the reference pulse for ATS-3 is provided by the sun. Given the orbital 
and the spin axis characteristics, there is an angle, /3, between the sun reference and the local 
vertical vector. This angle 0 and the time rate of change, are precomputed and stored on 
the ground in the processing equipment. It should be remembered that the Soumi camera 
on ATS-3 is constantly on and scans the earth disk and the darkness of space as well. The 
digital processing occurs on the ground. 





SPIN AXIS 


Figure 2. Camera scan geometry 

Again referring to figure 2, the camera .steps from 9° above the plane normal of the spin 
axis to 9° below the plane normal to the spin axis. For a typical elevation angle, 5, the 
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Figure 3. Misalignment effects on pictures, (a) Spin axis at orbit normal and no misalignments; (b) spin 
axis off orbit normal and no misalignments; (c) spin axis at orbit normal and bias misalignment; (d) spin 
axis at orbit normal and skew misalignment. 
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ground equipment initiates recording at a point which is a function of tlie precomputed 
angle, and the angular rate of change, p. During the sc,.n, the leading edge of the earth is 
observed as well as the entire earth disk. At the point where the trailing edge of the earth 
is obsared, the scan is terminated. The spacecraft then rotates, the camera is stepped the 
next discrete elevation angle, and the process is repeated. 

The result is a matrix consisting of 2400 scan lines (corresponding to the 1 8° field of view) 
and 8192 samples per scan line (the analog signal is sampled 8192 times over the duration 
of the recording). The recording is synchronized with the spin rate such that a 20** scan is 
obtained. The matrix is searched by a preprocessor program called EDGE, and EDGE 
extracts either the leading edge of the earth or the trailing edge of the earth, depending upon 
which is the most sharply defined. Prior to high noon, we process the right earth edge. 

After high noon, we process the left earth edge. (Since this is a visible system, only one edge 
per picture is processed.) 

The effects of spin axis orientation and misaligiunents are illustrated in figure 3 with four 
limiting cases. C!ase a depicts a purely circular synchronous orbit with the spin axis at 
orbit normal and no misalignments present between the principal and the geometric axes. 

In this particular case, we would see a perfectly centered earth throughout four points in 
the orbit, each separated by 90°. The effects of terminators are not illustrated in the figure. 
With a vis’ble system, terminators would be present within these four pictures. 

If no misalignments are present but the spin axis is perturbed off orbit normal (case b), 
the top of the earth would be cropped at <me point in the orbit. The bottom of the earth 
would be cropped one half orbit later. One quarter and three quarters of an orbit later, the 
earth would ' ; centered perfectly in the picture. 

In case c, an orbit-normal attitude is represented, and a bias misalignment, d, is present. It 
should be noted that if the total misalignment angle between the geometric axis of ^in and 
the principal axis of spin is a biased misalignment, then the camera plane, the X-principal 
axis, and the X-geometric axis are ail coplanar. With this particular geometry throughout 
the orbit, a constant cropping of eith-i the north or the south of the earth within the frame 
would occur. 

Finally, with case (d), an orbit normal attitude is represented, and the total misalignment angle 
is the skew misalignment an^e In this case, there is a stretching of the earth, a distortion, 
because the angle P is computed by assuming that there are no skew misalignments. If there 
is a skew misalignment, recording will be sooner than desired at the top of the picture and 
later than deared at the bottom of the picture. 

In general, the spin axis i not at orbit normal, and misalignments are present. The objective 
is to use the orbital information and tho» estimate the spin axis attitude angles as well as the 
misalignment angles ^ and 6. 
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The basic problem of estimation in PICATT can be outlined as follows* 


PICATT Estimation Problem 


Constants C 

n = 

Right ascension of ascending node 

CO 

= 

Argument of perigee 

M 

= 

Mean anomaly 

i 


Inclination 

a 

= 

Camera elevation angle 

€ 


Earth optical radius 

r 

= 

Spin period 


State variables X 


a, 6, 6 

Measurement Y 

03 = Sweep angle from initiation of recording to earth edge 

The standard orbital parameters are known constants. The semimajor axis and eccentricity 
are not included. PICATT actually does account for the complete orbital state in the 
mathematical model. At synchronous attitude, with a very low eccentricity, there is only 
second-order effect from those two parameters. We use the argument of perigee and the 
mean anomaly to describe the intrack position. In addition to intrack position, the inclina- 
tion is also required to define the orkntation of the local vertical relative to the ^in axis. 
The camera elevation angle is also assumed to be exactly known, as is the earth optical 
r .dius (the angular width of the earth as viewed from synchronous altitude, approximately 
17.4°) and the spin period (needed to convert the 8192 discrete digital counts to a sweep 
angle). These known parameters— the orbital information, camera elevation angle, earth 
optical radius, and spin pe::od— are all assumed to be known with zero uncertainty. 

A four-parameter state variable array includes the right ascension (a), declination (5), the 
skew misalignment angle (^), and the bias misaligrunent angle (d). The precession of die 
spacecraft due to solar radiation torque is ignored because picture information is obtained 
over only a quarter of an orbit. In other words, 6 hours of information, rather than 24 
hours, is obt-rined. Over the computational period in which the pictures are taken, it is 
assumed that a and 5 are inertially fixed. 

Finally, the measurement Y is the sweep angle from the initiation of the recording to the 
detection of either the leading edge or the trailing edge of the earth. This measurement is 
obtained by converting the edge counts to an angle of sweep, using the spin period, r. 

PICATT uses the weighted least-squares batch filter with memory and requires knowledge 
of the Bayes matrix, A^ (the inverse covariance matrix of error and the initial state estimate) 
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Vi 


= (aT wr' A. + Aj) ' (aT W-' y. + A-„' x.] 


where 


X. = 

1 

x-x 

1 

= 


Vi = 

Y-f(X,,C) 

^ = 

(3f/3X). 

W ‘ = 

Weighting matrix 

V = 

Bayes matrix 


The Xj is the residual between the unknown true state and a current best estimate of the 
state; x. is the residual between the initial state estimate and the current best estimate of 
the state; y. is the residual between the measured sweep angle Y and the predicted sweep 
angle, which is a function of our current best estimate of state as well as the known con* 
stants. The A matrix is the matrix of partial derivatives of the observation equation, and 
W is the weighting matrix. In PICATT, the data are weighted as a function of the computed 
standard deviation of the residuals between the measured and observed sweep angles. 

How well does PICATT work? Prior to the development of the PICATT program, the 
operational attitude determination program for ATS-3 was the ATBAY program, which 
used sun sensor information and polarization angle measurements (POLANG). The prob- 
lem with ATBAY was the POLANG. Real-time Faraday rotation measurements normally were 
not available. In addition, the stations tended to have biases that could not be estimated. 

It was felt that this was no better than a I ° accuracy system. 

Prior to the development of the PICATT program, we would perform an attitude maneuver 
on ATS-3 at approximately 2-month intervals. That was a very ineffleient procedure as far 
as performing attitude maneuvers-trying to erect the spin axis to orbit normal and then 
finding out that the spin axi** actually precessed in some other direction. As a result, we 
were maneuvering approximately every 2 months. With the advent of the PICATT program, 
the requirements for maneuvering ATS-3 have been reduced to a frequency on the order of 
one maneuver every 8 months. Because of the increased accuracy of PICATT, the efficiency 
has been significantly improved as far as the attitude control. 

The following is a comparison between the ATBAY program and PICATT, for the time per- 
iod covering August 1970 to December 1971. During this time, we were not receiving pic- 
tures on a weekly basis as we are now, so we did not have a large volume of PICATT solutions. 
Normally, the ATBAY solutions were computed out of phase. For this statistical study, 
we selected spin axis attitude estimates at points that are separated by no more than I or 2 
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days between ATBAY and PICATT. With this ground rule, there were 1 1 samples. The 
residuals are the difTerences between ATBAY and PICATT. It is assumed that PIC.ATT is 
correct. Listed below are the mean (ji) and the standard (5) deviations for these 1 1 samples; 

M 6 


Right ascension 

0.512 

13.486 

Declination 

-0.195 

0.225 

Pointing error 

0.496 

0.302 


The most important characteristic is the pointing error. With this statistic, the mean is 
almost 0.5° with a standard deviation of 0.3; 8°. The sum of the mean and standard 
deviation approach the 1° total error in the sun sensor/polarization angle method of 
attitude solution. 

This is a relative comparison; To find out the actual operational results on ATS-3, an inter- 
mediate computational program called the latitude scan program is used. This program takes 
the orbital information, the attitude and bias state estimate (a and 5 as well as the misalign- 
ment angles <ff and 0), and r^^d’ots the latitude that will be -ned in the center of the 
picture as a function of * lay. 

The first scale on figure 4 oa'nts right ascension, which ranges from 0° to 360°. Time 
of day can be related to rignt : cension in inertial space, since the earth is moving .'C>0° in 
one day. Likewise, !o d hi,"h n. jn as a function of time of year can be correlated as a 
function of right ascension. The solid line 'ndicates a prediction for observing the maximum 
north latitude scanned on the pictures as a function of time, using the PICATT state estima:e 
and the orbit parameters. Tne portion that is fairly flat indicates that the camera is scanning 
over the top of the earth. The sharply curved portions indicate that the camera is scaiming 
into the earth disk. The triangles represent latitudes whidi have been read off gridded 
weather pictures. It should be remembered that not only do we get the digital pictures, but 
we can also obtain a corresponding print by processing the video signals from the camera. 

It is seen that, the farther the spin axis is located from orbit normal, the more the scan cuts 
into the earth. At the point that the scan starts cutting into the earth, seeing lower and 
lower, and approaching the 55° north latitude limit, the better observability we have as far 
as the correctness of the PICATT solution is concerned. From the figure, it is seen that we 
are predicting very accurately where the cutoff occurs between the point where we scan over 
the top of the earth and the point where we start cutting into the center of the earth. Nor- 
mally, we find that right ascension is the most difficult parameter to define; we usually de- 
tect a bad solution by finding a horizontal shift in the prediction. We generally get fairly good 
agreement on the vertical scale or latitude. 

It has been indicated previously that crossHroupUng of attitude perturbations with the east- 
west stationkeeping maneuvers occurs. The last stationkeeping maneuver was performed 
December 5, 1973. Figure 5 presents the ATS-3 state since then. There is a general trend 
in right ascension and declination with some scatter. There is even more scatter in the total 
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Figure 4. Typical latitude scan plot. 

misalignment angle, as seen in figure 6 (total misalignment angle being the included angle 
between the geometric axis of spin and the principal axis of spin). The problem illustrated 
is traceable to the fact that we get only a quarter orbit of information, and the spin axis is 
fairly close to orbit normal. Being close to orbit normal, the PICATT program has difficulty 
relating the bias misalignment angle, 6 , to the declination state, S. There is a trade-off be- 
tween 9, which is the primary contributor toward the total misalignment angle, and the 
declination, 5. 

We currently are obtaining erratic results, this would improve if we were to have information 
over a full orbit. Of course, being a visible camera system, we cannot get that information 
over a full orbit. Included in figure 6 is a statistical study on the misalignment information. 
We are predicting a mean misalignment of 0.348° with a standard deviation of O.l 12° and 
a total misalignment angle on the order of 0.5°, although earlier in the development of the 
PICATT program, it appeared that the misalignment angle was on the order of 1°. Every 
time an east-west stationkeeping maneuver is performed, more fuel is used and the misalign- 
ments change. 

The PICATT program has been extremely successful, as far as maintaining the attitude on 
ATS-3 within the given requirements. I think that the program has more capability if it can 
be used on a satellite that has digital earth pictures over a complete orbit. This is possible 
with the Synchronous Meteorological Sateilite-A (SMS-A), launched May 16, 1974, which 
has an infrared capability and takes pictures over the complete orbit. We expect SMS to 
get much more consistent estimates of the misalignment an^es as well as the attitude angles, 
a and 5. 
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Figure 5. ATS-3 attitude time^iistory following December 5. 1973, 
stationkeeping maneuver. 
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Figure 6. ATS-3 Misalignment time-historv following Oecenrber 5, 
1973, stationkeeping maneuver. (Misalignment angle statistics: ft • 
0.348\ 0*0.112°.) 
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By intrack orbital error is meant a constant time adjustment, AT, that is applied to a set of 
efriiemeris data which is otherwise correct. The ephemeris data may be in the form of an 
orbit tape or in the form of orbital elements with an associated orbit generator. The AT is 
simply added to the time before the ephemeris routine is accessed. It is imphcit here that 
AT is a constant throughoui the pass of data that we are considering, where the pass of data 
is typically a fraction of one orbit. 

In terms of Keplerian orbital elements, we are applying an adjustment to one of the six orbi- 
tal elements, and the other five are assumed to be correct. Why would we want to assume 
that flve of the ax orbital elements are correct? Those who are familiar with orbit determin- 
ation problems understand that there are cases, particularly with a predicted orbit, where the 
predominant source of error will be an intrack error. We have seen numerical examples from 
real, predicted orbit tapes compared with later definitive and more accurate orbit tapes, which 
show that as much as 99.9 percent of the orbit error in a predicted tape can be removed, 
simply by applying a constant time adjustment throughout an orbit. Table I shows an ex- 
ample of intrack orbit error for the Small Scientific SatelUte-A (SSS-A), which has an apogee 
height of 26,500 km, a perigee height of 220 km, and a period of 7 hours and 20 minutes. 

The enor in the predicted orbit tape is determined by comparison with the definitive tape. 
The predicted tape is accessed at a time about 2 weeks beyond the available data used in 
the predicted tape. 

Table 1 

Example of Intrack Orbit Error 



Error in Predicted 

Optimum Time 

Error After Time 


Tape (km) 

Adjustment (s) 

Adjustment (km) 

Near Apogee 

123. 

59.13 

1.83 

Near Perigee 

586. 

59.24 

0.50 


The primary motivation for determining this intrack adjustment is to improve the accuracy 
of our attitude determinations, oarticularly in cases where we are forced to determine an 
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attitude in near real time within an hour of the time the data are received, in tnese cases, 
we have to use a predicted orbit tape; we do not have time to wait for a definitive orbit tape 
to be generated. Potentially, this technique has the capability of improving orbit determina- 
tion for other users as well or of attaining the same orbit accuracy that we have now, but 
using less orbit data. That has not yet been done, but we are working on combining the 
orbit and attitude problems, that is, processing orbit tracking data with earth and sun sensor 
data in one system and thereby improving both t/.e orbit accuracy and the attitude accuracy 
with the available data. 

Figure 1 explains these attitude sensors, which have been used on at least four different 
missions-the Radio Astronomy Explorer-2 (RAE-2), the Interplanetary Monitoring Platfonn 
(IMP), the Small Scientific Satellite (S^ ), and the Atmospheric Explorer (AE)-and which 
are planned to be used on the Synchronous Meteorological Satellite (SMS) and the Communi- 
cations Technology Satellite (CTS). What all these missions have in common is an earth sensor 
telescope of some type, mounted at an angle to the spin axis so that the earth sensor scans a 
cone; if this cone intersects the earth, then the earth sensor will be triggered. The telescope 
may be sensitive to either infrared or visible light. In addition, there is a sun sensor on the 
spacecraft with a slit parallel to the spin axis; when the plane of that slit crosses the sun, the 
sensor triggers and also measures the angle between the spin axis and the sun direction. 

The raw telemetry includes the angle between the spin axis and the sun, (/3); the time that the 
sun sensor slit plane crossed the sun; the times that the earth sensor triggered on and off; and 
the inertial spin period, as defined by the time between two successive sun sightings. 

As shown by figure 2, it is easier to visualize the information if it is considered in terms of 
the geometric parameters which it defines. It happens that all the information in a single 
frame of data defines only three angles in space. One of them, of course, is the sun angle, 
the angle between the spin axis and the sun, because it is measured directly. The other two 
geometric parameters are the dihedral angles labeled A^^ and A^^^^ in figure 2: Aj^, is the 
dihedral angle from the plane of the spin axis and the sun to the plane defined by the spin 
axis and the horizon vector, that is, the vector from the spacecraft to the horizon at the earth- 
in triggering; A^ is the same thing for the earth-out triggering. It is worth noting that these 
two horizon vectors are unknown quantities. Even if the vector from the spacecraft to the 
earth is known, the vector from the spacecraft to the horizon crossing points would not be 
known. The data define two dihedral angles measured with respect to unknown vectors. 

We will now examine the resolution of these devices for the missions considered. The data 
are of course digitized before we receive them, so the value of the least significant bit is a lower 
limit on the resolution of the sensor. This is not to be confused with the accuracy of the 
sensor, because there may be systematic errors much larger than the resolution. 

For the sun angle, the least significant bit typically has a value ranging from 0.25° to 1°, 
making the sun sensor a relatively coarse sensor. The earth-in and earth-out measurements 
are somewhat more sensitive. For these, the resolution depends on the clock rate of the 
spacecraft in relation to the spin period, and the resolution ranges from 0.01° to 0.7°. We 
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are mainly interested in missions where the resolution is closer to 0.01° because, in that 
case, the earth sensor is potentially a rather sensitive device, assuming that we can remove 
systematic errors, which may be as large as 1° or more, and which will be discussed later. 

This is what the data lock like in a single frame, and now 1 want to define the unknowns in 
this problem. The first is the attitude of the spacecraft, because that is what the sensor 
was put on board to determine. The attitude can be described by a two-element state vector, 
right ascension (a) and declination (5), if we assume that the attitude is constant and there 
is no nutation. I am assuming throughout this presentation that any nutation in the space- 
craft is negligible. Therefore, attitude is a two-element state vector. The second unknown 
is the intrack time adjustment, At, which is the primary topic of this paper. We are assuming 
that we have a source of orbital ephemeris, which is correct, with the possible exception of 
this intrack time adjustment. 

So these are three primary state oarameters: two for attitude and one for time adjustment. 

In the ideal case, those would be the only three unknowns in the problem, and the problem 
would be relatively simple. In practice, it has been found that, on all of the missions we have 
supported, there are significant systematic biases in the sensors that have to be removed 
when the data are processed in order to meet the attitude requirements of the mission, which 
may be, for example. ±1° for attitude. 

Therefore, some additional parameters must be determined from the data. In general, the 
unkiiowns include the elevation of the earth sensor with respect to the spin axis, the aximuth 
of the earth sensor with respect to the sun sensor, and the elevation of the sun sensor with 
respect to the spin axis, that is, a bias in the measured sun angle. In addition, there is a 
possible earth sensor triggering threshold or sensitivity error. These earth sensors do not have 
a very narrow field of view; the field of view may be as wide as 3° in diameter. If the sensor 
threshold is not accurately known, there may be an uncertainty of several degrees as to where 
the sensor is pointing at the time it triggers. Finally, there is the possibility of a constant time 
delay on either the earth-in or the earth-out triggering due to electronic delay between the 
time the event occurs and the time it is recorded. 

In principle, all of these quantities can be measured on the ground before the spacecraft is 
launched. In practice, they are subject to change. The alignments, of course, could change 
due to thermal distortion of the spacecraft. It is even more likely that the apparent align- 
ments with respect to the s-'in axis would change because the spin axis shifts with respect 
to the geometric body axis, for example, due to uneven fuel usage between fuel tanks. Also, 
the electronic parameters can change, for example, if the temperature of the electronic com- 
ponents changes. 

To model all these sources of error, it is necessary' to introduce five additional angular 
parameters. These are biases with nominal values of zero. There is a bias on the earth sensor 
elevation and a bias on the earth-in or the earth-out rotation angle. This includes both the 
effect of an azimuth offset between the earth and sun sensors and a possible difference be- 
tween the time delays for the earth-in and the earth-out triggerings. There is also a bias on 
sun angle and a bias on the apparent angular radius of the earth as seen from the spacecraft. 
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The bias on the angular radius of the earth is intended to correct for the earth sensor trigger- 
ing threshold. Without discussing the details, it can be shown that, for an earth sensor with 
a circular field of view, any constant triggering threshold can be exactly compensated for by 
adding or subtracting a constant bias on the apparent radius of the earth. 

Thus, the total number of unknowns includes five angular biases and three primary state 
parameters, for a total of eight parameters that have to be determined. However, a frame of 
data includes only three observables: sun angle and eartn-in and earth-out angles. Clearly, 
on the basis of a single frame of data, we could determine, at most, three of these unknowns. 
To have any hope of determining ali eight unknowns, we need more than one frame of data, 
and the frames have to be independent in some sense. This is where the real problem occurs. 

Typically, we do have a large number of frames of data, but the sun angle may be constant 
throughout the entire block because, as mentioned previously, le sun sensor is a relatively 
coarse sensor, and the sun angle is changing very slowly. So it is not uncommon for the 
measured sun angle to be constant throughout the pass, in which case, regardless of the num- 
ber of frames of data, there is only one actual observable for the sun angle. 

The earth-in and earth-out angles are more useful, because they do vary with the spacecraft 
position. Still, two frames of data taken at nearby positions in the orbit will be redundant. 
Speaking qualitatively, in order to determine all eight of these parameters, it is clear that 
we need a significant fraction of an orbit of data in order to have independent observables 
and not just the three observables that occur in one frame. 

We have developed a program called OABIAS, which processes data of the type described and 
determines a state vector, including the unknowns listed below: 

• Attitude (a, 6) 

• Earth sensor elevation bias 

• Bias on earth-in angle, Ajj, 

• bias on earth-out angle, A^ 

• Bias on angular radius of earth 

• Bias on sun angle, 0 

• In-track orbit time adjustment. At 

The program is a standard, weighted least-squares recursive estimator. Any number of the 
above fisted parameters can be fixed at constant values and not determined. 

The program works as a standard recursive estimator: the state vector is used to predict the 
observables, the residual is computed for each observable, and the partial derivatives of each 
observable with respect to each element in the state are computed. Then the residuals and 
partial derivatives are used to update the state vector. The partial derivatives can be computed 
analytically for every case except the case of interest here, the intrack orbital time adjustment. 
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If the partial derivative of any arbitrary observable is considered with respect to the time 
adjustment, At, it is equal to the partial derivative of the observation with respect to the 
spacecraft position vector, multiplied by the partial derivative of the spacecraft position 
with respect to At: 

9 (observation) 9 (observation) 9R 

9At 9R 9At 

9 (observation) 

= • V 

9R 

The derivative of the spacecraft position with respect to At is one that cannot be computed 
analytically, because we do not have an analytical expression for spacecraft position as a 
function of time. However, it is not necessary, because we can get the velocity (V) cf the 
spacecraft from the orbit tape. 

In practice then, we analytically compute the derivation of each observation with respect to 
then multiply that by the velocity vector obtained from an orbit tape. The important 
point here is that this method is not restricted to any particular type of orbit. We are not 
assuming, for example, a Keplerian orbit. Any orbit that can be described by an orbit tape 
can be handled correctly using this technique. 

Figure 3 shows some results obtained using simulated data for the Communications Tech- 
nology Satellite, which is scheduled to be launched next year. The data were simulated for 
transfer orbit, with a perigee of 190 km and an apogee of 36,000 km. The attitude in this 
case is pointing 49° below the plane of the orbit; the earth sensor is an IR sensor mounted 
at 85° from the spin axis. These facts together imply that the earth sensor will scan the 
earth only during the indicated portions of the orbit. For the rest of the orbit, the earth 
sensors will miss the earth; there will be no useful data from the earth sensor for those 
periods. 

Both of these sections of the orbit were simulated and the data were combined. The data 
includes 100 frames evenly spaced over 240 minutes over both of those segments of the 
orbit-40 minutes near perigee and 200 minutes near apogee. Gaussian noise of 0.012° is 
applied to the earth rotation angles, the earth-in and earth-out angles. That corresponds to 
the clock rate expected for the spacecraft. 

We were attempting to determine the complete eight-element state vector, so biases of 1° 
were applied to each of the angular parameters that describe the bias on the sun angle, the 
earth-in and earth-out azimuth angles, the earth sensor elevation, and the apparent angular 
radius of the earth. For the time adjustment, an error of 60 seconds was applied to the 
ephemeris data, which corresponds to a very large, 6° error in spacecraft position at perigee, 
which would make the perigee data virtually useless for attitude determination without 
correcting for i» But the time adjustment corresponds to an error of only 0. 1 5° true anomaly 
at apogee. 
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Figure 3. Example of position determination 
with simulated data for the Communica- 
tions Technology Satellite. 

The results for these simulated data are shown in table 2. The biases are initially cfimated 
at zero, because they are unknown quantities; the initial attitude was obtained fr^m a deter- 
ministic processing of the data, before correcting for biases. It can be seen that this initial 
attitude is more than 2 ^ off, which would violate mission constraints for this mission, since 
we have a 1° attitude accuracy requirement. The final result f>om the program shows that 
all of the unknowns are determined to an accuracy of 0.05°, and the intrack time adjust- 
ment is determined to an accuracy of one-half second. This indicates that, first of all, the 
program works. Secondly, it means that the problem is feasible. That is, we really should be 
able to determine all eight of these state parameters, assuming that we have a sufficient 
amount of data and that there are no substantial systematic errors that have not been con- 
sidered in the state vector. 

This brings us to the case of the real data. We have processed real data from four spacecraft— 
IMP, S’ , RAE-2, and AE— but the results cannot be presented in this form because the true 
state is not known for any of these spacecraft. In fact, since the attitude changes, or is 
changed, from one day to the next, we do not generally have more than one pass of data to 
examine to define a given state vector. 
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Tabk 2 

OABIAS Results* 


Parameter 

True 

State 

— 

Initial 

Estimate 

Final 

Result 

1 

330.00 

329.50 

330.02 

.Attitude 1 g 

-21.75 

-23.88 

-21.80 

Earth Sensor Elevation 




Bias 

-1.00 

0.0 

-1.01 

Earth-in Azimuth Bias 

-1.00 

0.0 

-0.99 

Earth-out Azimuth Bias 

-1.00 

0.0 

-0.9b 

Bias on Angular Radius 




of Earth 

1.00 

0.0 

1.01 

Sun Angle Bias 

-1.00 

0.0 

-0.99 

Intrack Time Adjustment 

-60.00 

0.0 

-59.59 


* All values arc in degrees, except intrack time adjustment, which is in seconds. 


Thus, only three things can be said about the real data: either it can or cannot be fitted, or 
we can find an infinite number of state vectors which would fit the data equally well to 
within t:io noise level. It happens that this third cause is not uncommon for the missions 
that we examined. The reason is that, when trying to determine an eight-element state vec- 
tor on a small section of an orbit, it is to be expected that tlie entire state vector will not be 
determinable. 

However, we have generally found that the results are consistent with what we would expect, 
based on simulation. That is, where there are enough data to find the state vector, it can be 
determined uniquely. Where there is not enough information, we can determine any number 
of state vectors that would fit the data. We have not, as a rule, encountered the other prob- 
lem, that of not being able to find any state vector that fits the data; this would '. ..licate the 
presence of some systematic error that we have not modeled. 

It should be pointed out that, if the misalignment parameters and biases could be eliminated 
and the state vector thereby reduced from eight elements to three elements, it would ob- 
viously be a simpler problem. It would take much less data to determine the state. In fact, 
since there are three observables in a single frame of data, and as there would be only three 
'nknowns in that case, the time adjustment and the attitude could be determined, based on 
just a single frame of data. But. based on our experience, it is not feasible to ignore all of 
these sys'ematic bias parameters. 
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DISCUSSION 

VOICE: Is this technique being used operationally? 

SHEAR- It has been applied operationally to all four spacecraft mentioned: currently, the 
operational work is concentrated on AE-C. It has improved the attitude determination 
accuracy for those missions. It is not used to process every pass of data. It is used on the 
initial passes of data in the mission to try to determine the biases that will be used througli* 
out the mission. 

VOICE: Are the simulation results that you've shown for the same static attitude estimator? 
When you generated data, did you use the model you have in your estimator? 

SHEAR: That's right. It's the same model. There are no systematic errors tliat aren't 
accounted for. 

VO/Ct. There are no attitude dynamics? 

SHEAR: Correct, there are no attitude dynamics. 

VOICE: In the case of the real data, the accuracy of your intrack time adjustment could 
be checked against orbit tracking data, couldn’t it? 

SHEAR: It is possible to get an independent confirmation of the intrack orbit error. I 
think our problem has been mote otten that we can’t determine all eight state parameters 
including the time adjustment; consequently, if the time adjustment does not agree with the 
orbit data, we don't know whether we've got the right answer for one of those biases or not. 

I thinK that it is feasible to determine the intrack time adjustment on the spacecraft. The 
has an orbit that is fairly similar to the orbit 1 just showed you for the simulated example. 
In that case, I think it’s feasible, but we haven’t processed data extensively yet. 

VOICE: 1 have worked on some of the problems of bias determination and tiic biases are 
usually so highly correlated with what you’re trying to measure that they can't be separated. 
Is your simulation realistic? 

SHEAR: Well, there are different geometry cases to be considered. I could point out. for 
example, a circular orbit. For any circular orbit, we cannot determine this complete eight- 
element state vector, regardless of the attitude of the spacecraft, the sensor mounting angles, 
or the amount of data. It can be proven that there is a perfect correlation between the in- 
track error and the sensor misalignments; we can find an infinite number of state vectors 
that will fit the data. 

Referring back to figure 3, if we were to delete either one of these pas.ses of data. I don't 
think it would be feasible to do the problem; the complete eight-element state vector could 
not be determined with either pass alone. Of course, if we could reduce the stale vector 
to fewer elements, it could be done. But we were not able to solve for all eiglit elements, 
even on simulated data, with just one of those sections. 
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VOICE: As a practical problem, can you get track-in data when you're at perigee? 

SHEAR: Probably not directly at perigee, for such a low perigee. But it should be possible 
to get most of these data. The spacecraft, which has the same orbit except for a sliglitly 
lower altitude, frequently obtains a pass of data which covers most of what is indicated here 
as the perigee pass. 

VOICE: For your simulated data, where did the sun lie with respect to the orbit flight? 

SHEAR: The sun is shown in figure 3. but i didn't mention it. The sun is 32“^ above the 
orbit plane, and it’s rouglily in the direction shown in the figure. 



